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I.   INTRODUCTION 

A.   HISTORY 

The  existence  of  clouds  and  fog  in  many  regions  of  the 
earth  presents  a  formidable  problem  to  the  designer  of  an 
optical  communication  system  whose  transmission  channel  is 
the  atmosphere .   A  scattering  medium  can  inhibit  system  per- 
formance by  inducing  beam  spread,  dispersion  in  angle-of- 
arrival ,  degradation  of  spatial  coherence  and  dispersion  in 
time  and  frequency  of  the  modulated  optical  beam.   The  develop- 
ment of  an  optical  communication  system  for  atmospheric 
applications  requires  an  accurate  knowledge  of  the  effects  of 
scattering  on  light  propagation.   Numerous  studies  have  been 
made  in  an  attempt  to  provide  this  knowledge.   Different 
aspects  of  the  problem  have  been  assaulted  using  various 
mathematical  models.   Monte  Carlo  computer  simulation  has 
been  used  by  Bucher  [1],  Plass  and  Kattawar  [2]  through  [5], 
Junge  [6]  and  Danielson,  Moore  and  van  de  Hulst  [7],  to  mention 
only  the  few  used  extensively  in  this  work.   Equally  prominent 
attempts  have  been  made  using  analytical  methods  such  as  those 
by  Arnush  [8],  Gordon  [9],  Stotts  [10],  Hansen  [11],  Ishimaru 
[12],  Kennedy  [13]  and  Lutomirski  and  Yura  [14-].   Complicated 
numerical  techniques  have  been  used  only  by  Dell-Imagine  [15] 
and  Zachor  [16].   Each  attempt  has  contributed  to  the  stock- 
pile of  knowledge  required  but  additional  problems  remain  to 
be  investigated.   Many  books  have  been  published  on  the  subject 
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of  light  scattering.   Those  found  most  useful  to  this  work 
were  van  de  Hulst  [18],  Kerker  [19],  Chandrasekhar  [20]  and 
Deirmendjian  [21].   Interested  parties  and  reasons  for  their 
interest  are  briefly  discussed  in  the  following  paragraphs . 

B.   PURPOSE 

Current  navy  operational  communications  systems  suffer 
from  a  number  of  problems .  There  is  no  operational  system 
which  is  not  subject  to  jamming,  intercept,  spoofing  and 
direction  finding.  In  addition,  it  appears  that  optical 
communications  systems  have  great  promise  in  solving  these 
problems  for  many  applications  [17].  Any  information  that 
could  be  used  in  evaluating  such  systems  is  desired. 

The  Navy's  effort  to  create  a  worldwide  satellite-to-sub- 
marine communication  network  is  another  rapidly  progressing 
area  requiring  information  of  the  nature  being  investigated 
in  this  thesis.   In  this  situation,  light  propagation  through 
water  further  complicates  the  matter. 

Among  other  parties  that  may  require  information  on  light 
propagation  are  the  United  States  Coast  Guard  in  their  aids 
to  navigation  system,  NASA  in  satellite  communications  and 
researchers  trying  to  determine  the  composition  of  the  atmos- 
phere by  analyzing  scattered  radiation. 

The  purpose  of  this  work  is  to  characterize  both  the 
spatial  and  temporal  aspects  of  a  light  beam  propagating 
through  a  scattering/absorbing  medium  using  an  analytically 
verified  Monte  Carlo  model  of  the  system.   Models  were  created 
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and  tested  against  each  other  until  an  adequate  level  of 
confidence  in  each  was  obtained.   The  models  were  then  used 
to  study  more  specific  questions  concerning  time  and  spatial 
spread.   More  explicitly,  the  models  were  used  to  observe 
what  effect  different  phase  functions  had  on  the  temporal  and 
spatial  aspects  of  light  scattering.   It  was  believed  that 
the  pronounced  forward  peak  of  some  phase  functions  would 
generate  different  time  and  spatial  character  than  a  moder- 
ately peaked  phase  function.   An  investigation  attempting  to 
verify  this  belief  was  an  important  part  of  this  work. 

Spatial  characterization  of  light  passing  through  a  cloud 
of  finite  thickness,  and  transition  from  the  forward  scatter 
to  diffusion  region  was  also  investigated. 

Once  confidence  was  established  in  the  models  used,  the 
general  areas  of  agreement  or  disagreement  were  studied  in 
an  attempt  to  specify  when  one  model  may  be  used  more  effici- 
ently than  another  for  generating  information  concerning  a 
specific  situation.   The  drawbacks  and  strong  points  of  each 
model  type  have  been  determined  and  mentioned  where  appropriate 
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II.   THEORETICAL  PRELIMINARIES 

A.   INTRODUCTION 

The  purpose  of  this  section  is  to  review  the  various 
theoretical  formulations  currently  used  to  characterize 
optical  propagation  in  a  single  and  multiple  scattering  medium. 
A  complete  coverage  of  any  one  of  the  topics  in  the  following 
sub-sections  is  most  certainly  not  contained  here  but  suffici- 
ent references  are  provided  for  the  reader  inclined  to 
investigate  any  topic  in  more  detail.   Only  those  details 
related  closely  to  this  work  are  described  in  any  detail  and 
even  then  only  the  important  conclusions  are  mentioned  in 
many  cases. 

The  order  of  presentation  is  as  follows .   The  principal 
characteristics  of  the  scattering/absorbing  medium  are  men- 
tioned briefly,  thereby  defining  frequently  used  terms.   A 
brief  summary  of  Mie  theory  is  included  to  present  the  valuable 
insight  necessary  in  investigating  the  problems  posed. 
Radiative  transfer  theory  is  mentioned  in  passing  because  the 
solution  of  the  complicated  radiative  transfer  equation  is 
essentially  the  subject  of  many  analytical  attempts  at  solving 
the  scattering  problem.   This  section  then  probes  into  mathe- 
matically modeling  the  geometry  of  a  real  atmosphere .   Three 
general  categories  of  models  are  discussed. 

Analytical  models,  which  made  various  assumptions  allowing 
closed  form  expressions  for  certain  characteristics  of  the 
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medium,  are  discussed.   The  approximations  in  such  models 
normally   detract  from  the  total  worthiness  of  such  an 
approach. 

The  discussion  then  progresses  into  the  use  of  numerical 
methods  in  evaluating  analytical  solutions  as  well  as  in 
evaluating  the  radiative  transfer  equation  directly.   This 
often  overlooked  powerful  technique  has  been  used  by  theore- 
ticians because  of  its  ease  in  solving  complicated  integro- 
differential  equations  for  which  exact  solutions  have  been 
proven  impossible. 

Monte  Carlo  methods,  used  predominantly  in  this  work,  are 
then  discussed.   This  technique  is  used  extensively  due  to  its 
ease  in  adaptation  to  odd  geometries  of  the  scattering  medium. 
In  this  work  this  benefit  was  used  only  to  model  finite 
clouds  in  an  otherwise  homogeneous  atmosphere. 

This  section  is  completed  with  a  general  discussion  of 
the  results  predicted  by  theoretical  insight  into  the  problem 
of  light  transfer  in  a  random  medium. 

B.   PRINCIPAL  CHARACTERISTICS  OF  SCATTERING/ABSORBING  MEDIA 

This  section  divides  the  characteristics  of  a  scattering 
medium  into  three  sub-sections.   The  first  defines  the  para- 
meters of  the  atmosphere  which  during  the  course  of  this  work 
were  considered  constants.   The  second  considers  the  types  of 
scattering  encountered  in  a  scattering  medium  and  the  third 
briefly  describes  the  significance  of  the  phase  function  in 
scattering  theory. 
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1 .   Parameters 

In  general,  a  medium  may  exhibit  both  scattering  and 
absorption.   Some  of  the  energy  of  a  collimated  incident  beam 
on  a  particle  will  be  scattered  over  all  possible  directions; 
the  rest  will  be  absorbed  by  the  particle  and  lost  from  the 
radiation  field.   In  the  single  scattering  problem  the  inten- 
sity of  radiation  at  any  point  along  the  direction  of  plane  wave 
propagation  is  given  by: 

-(K    )R         -(K  .  +K    )R 
I(R)  =  I(0)e     exr    =  I(0)e     a£)S   sca   ,  (1) 

where  the  absorption  loss  and  scattering  loss  together  is 
called  extinction.   The  ratio   of  scattering  coefficient  to 
the  extinction  coefficient  is  called  the  albedo  of  single 
scatter.   Thus, 


%    =    K^K     '  °  <  %  <  1  (2) 

abs   sca 


It  is  sometimes  more  convenient  to  discuss  optical  dimensions 
in  terms  of  mean  free  paths  of  photons,  i.e.,  the  average 
distance  between  collisions.   The  mean  free  path  is  often 
called  an  extinction  length  and  is  given  by  the  reciprocal  of 
the  extinction  coefficient.   Optical  thickness  of  a  given 
medium  is  a  dimensionless  number  representing  its  thickness 

in  extinction  lengths  or  mean  free  paths. 

The  transmittance  of  the  atmosphere  is  the  fractional 

intensity  remaining  in  a  beam  after  traversing   a  path  R  units 
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long.   Explicitly,  with  K     in  km"   and  R  in  kilometers, 


ICR)  _  „  "KextR 


(3) 


iToT 

The  term  visibility  is  defined  as  [22] 


=  K— I  ln  702  =  "K—  '  (4) 

ext  ext 


thus  the  transmittance  for  a  path  length  just  equal  to  the 
visibility  is  two  per  cent.   In  many  atmospheric  light  propaga- 
tion problems  the  albedo  is  close  to  unity  and  K     in 

J  ext 

the  above  equations  can  be  expressed  by  K    with  little  error. 
H  *  J      sea 

In  a  real  atmosphere  there  is  a  contribution  to  scat- 
tering by  small  particles  as  well  as  a  contribution  by  large 

particles.   K     then  is  defined  in  even  more  detail.   The 
r  sea 

following  paragraphs  present  some  of  these  details. 
2 .   Types  of  Scattering 

Atmospheric  light  scattering  can  be  classified  in  two 
general  categories.   Small  particle  scattering,  where  the 
particle  radius  is  much  scalier  than  the  wavelength  of  the 
incident  light,  is  called  Rayleigh  scattering,  and  other 
scattering  from  larger  particles  is  termed  Mie  scattering. 
The  terminology  of  this  work  refers  to  KR    as  the  portion  of 
the  scattering  coefficient  due  to  Rayleigh  scattering  and 
}L..      as  the  portion  due  to  Mie  scattering.   Therefore, 
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K^  =  K,+K„.+Kr),  (5) 

ext    abs    Mie    Ray' 


Mie  scattering  is  often  called  particulate  scattering  while 
Rayleigh  is  often  referred  to  as  molecular  scattering.   A 
ratio  of  particulate  to  total  scattering  is  defined  for  use 
in  this  thesis  as: 

K  . 
RPT  =    ^e _  (6) 

KMie  KRay 

Scattering  may  be  isotropic  (scatters  in  all  direc- 
tions equally)  or  anisotropic  (scatters  as  a  function  of 
angle  off  incident  direction).  The  latter  only  is  used  here. 
The  problem  of  defining  scatter  for  single  particles  is  well 
understood.  The  theory  of  scattering  has  been  extended  to  a 
medium  of  many  equal  sized  particles,  a  monodispersion,  and 
to  a  medium  of  many  different  sized  particles,  a  polydisper- 
sion  [18,  21]. 

When  a  photon  undergoes  only  one  collision  in  tra- 
versing a  scattering  medium,  single  scattering  has  taken  place. 
On  the  other  hand,  when  a  photon  undergoes  more  than  one 
collision,  multiple  scattering  has  occurred.   This  work  in- 
volves independent  single  and  multiple  scattering  in  spherical 
polydispersions ,  where  no  scattering  event  affects  other 
scattering  events.   For  further  details  related  to  scattering 
parameters,  types  of  scattering  and  other  more  specific  topics 
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the  reader  is  referred  to  an  excellent  text  on  the  subject, 
McCartney  [22].   The  following  paragraphs  focus  on  the  angular 
dependence  of  scattering  for  scattering  of  different  types . 
3 .   Phase  Function 

The  phase  function  expresses  in  a  formal  manner  the 
angular  dependence  of  scattering.   The  phase  function,  denoted 
here  by  P(0),  is  defined  by  van  de  Hulst  [18]  as  the  ratio  of 
the  energy  scattered  per  unit  solid  angle  in  a  given  direction 
to  the  average  energy  scattered  per  unit  solid  angle  in  all 
directions.   This  definition  requires  that  the  integral  of 
the  phase  function  be  normalized  to  unity,  which  is  to  say  that, 


Jj       P(9)  Sin 


1~  /  /   P(0  )  sin  0  d0d0  =  1  (7  ) 


When  expressed  in  this  way  the  phase  function  is  a  scalar 
function.   In  a  more  complex  analysis  of  scattering  theory, 
there  are  phase  functions  of  different  types  for  different 
polarizations  of  light.   These  functions  form  elements  of  the 
normalized  scattering  matrix  which  is  used  to  transform 
Stoke' s  parameters.   Detailed  discussion  of  the  parameters  is 
given  by  Refs.  18  and  21.   In  Appendix  B  a  program  to  evaluate 
the  Stoke' s  parameters  for  the  spherical  polydispersions  used 
in  this  work  is  presented  and  explained.   The  scalar  phase 
function  is  obtained  by  averaging  the  two  Stoke 's  parameters 
obtained  using  this  computer  adapted  Mie  theory. 

In  many  cases  the  phase  function  is  easily  approximated 
as  a  closed  form  function  of  scatter  angle.   In  these  cases, 
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the  Henyey-Greenstein  function  is  often  used  as  the  approxi- 
mating function  because  of  its  simple  form.  For  the  Henyey- 
Greenstein  phase  function  (Figure  1) , 


P(9)  =  ^~G  ) t-t-  ,  -1<G<1        (8) 

4ir(l+G  -2Gcos0) 


where  G  is  a  parameter  which  equals  <cos6>  for  the  phase  func- 
tion.  A  G  value  between  .80  and  .87  approximates  real 
atmospheric  phase  functions  quite  well  in  most  cases.   However, 
the  following  phase  functions  were  not  represented  well  using 
the  Henyey-Greenstein  phase  function  no  matter  what  parameter 
was  used. 

Phase  functions  used  in  this  work  were  of  Water  Cloud 
C.2  [21]  at  a  wavelength  of  .53  microns  (Figure  2)  which  was 
calculated  using  the  program  of  Appendix  B,  and  NOSC  Fog 
(Figure  3)  which  was  calculated  in  a  similar  manner  using  an 
experimentally  determined  particle  distribution  near  San  Diego, 
California.   For  further  discussion  of  phase  functions,  see 
"Methods  of  Simulation"  later  in  this  thesis.   References  18 
and  21  study  the  phase  function  and  its  generation  extensively. 

In  general  the  phase  function   is  quite  symmetric  in 
the  forward  and  backward  direction  for  particles  small  compared 
to  wavelength.   For  large  particles  the  phase  function  has  an 
increasingly  complicated  dependence  on  angle  of  scatter. 

As  one  might  imagine,  the  scalar  phase  function  plays 
a  significant  role  in  single  and  multiple  scattering  problems. 
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Henyey-Greenstein   Phase    Function 
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Water  Cloud  C.2  Phase  Function  [21] 
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NOSC  Fog  Phase  Function 
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Much  attention  has  been  focussed  on  defining  the  conditions 
under  which  the  entire  Mie  series  can  be  approximated  by 
analytic  expressions  such  as  the  Henyey-Greenstein  function 
mentioned  previously  [1,  7,  11,  18].   This  thesis  is  part  of 
the  quest  to  define  these  conditions. 

C.   MIE  SCATTERING 

Mie  theory  considers  the  scattering  and  absorption  which 
occurs  when  an  electromagnetic  wave  is  incident  on  a  spherical 
particle.   It  generates  desired  parameters  of  the  atmosphere, 
namely  the  extinction  and  scattering  cross  sections,  as  well 
as  the  phase  function.   Mie  theory  is  extremely  versatile  in 
that  it  can  be  applied  at  any  particle  size  to  wavelength 
ratio  and  can  be  extended  to  many  particles  and  more  important- 
ly polydispersions .   Appendix  B  explains  the  adaptation  of  Mie 
theory  to  machine  computation  for  collecting  scattering  data 
for  polydispersions  of  a  given  particle  distribution.   As 
mentioned  before,  Mie  theory  was  used  to  generate  optical 
parameters  of  the  Water  Cloud  C.2  [21]  model  at  a  wavelength 
of  .53  microns.   It  is  prudent  at  this  time  to  define  the  Mie 
scattering  parameter  as 


x  =  ^r  (9) 


where  r  is  the  radius  of  the  particle  and  X    is  the  wavelength. 

For  the  reader  interested  in  complete  derivations  of 
Mie  theory,  see  Chandrasekhar  [20],  Sekera  [23],  Born  and  Wolf 
[2H]  or  van  de  Hulst  [18]. 
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D.   RADIATION  TRANSFER 

Any  discussion  of  light  propagation  through  the  atmosphere 
would  not  be  complete  without  at  least  a  mention  of  the 
radiative  transfer  equation.   This  equation  models  mathemati- 
cally what  is  happening  as  light  transverses  a  scattering 
medium.   The  radiative  transfer  equation  in  its  simplest  form 
is  [10], 


(y  ^  +  1)  I  (t,  r,    u,  0,  t)  =  ^oIo  (T,  ?,  y,  0,  t) 


o 


2tt  1  (10) 

r f  PC6,  0;  6',  0')  I  (t,  r,  y',  0',  t)  dy'd0' 


where  I  is  the  scattered  radiance,  I   is  the  radiance  due  to 

'   o 

the  distributed  source  produced  by  the  incident  beam  transvers- 

ing  the  region,  y  =  cos8 ,  t  =  K   .  z  is  the  optical  thickness, 

r  are  the  spatial  coordinates,  0,  0'  are  angular  coordinates, 

&>   is  the  albedo,  P  is  the  scalar  phase  function  and  t  is  time. 
o  r 

Approximations  can  be  made  to  solve  this  equation  for 
I(t,  r,  y,  0,  t)  in  closed  form.   The  equation  has  basically 
two  components.   One  represents  the  field  caused  by  the 
incident  spatially  distributed  field  and  the  other  represents 
the  field  scattered  out  of  the  direct  beam  but  redirected  back 
into  the  same  direction  by  other  scattering  events.   For  a 
complete  derivation  and  explanation  of  this  equation, 
Chandrasekhar  [20]  is  suggested.   Under  specific  conditions, 
approximate  solutions  to  this  equation  can  be  found.   The 
next  section  mentions  a  few  of  the  efforts  made  toward  finding 
useful  solutions. 
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E.   ANALYTICAL  METHODS 

As  mentioned  previously,  several  authors  have  attempted 
to  provide  knowledge  of  the  effects  of  particulate  multiple 
scattering  on  light  propagation  by  applying  analytical  methods 
[8-13,  20,  25,  26].   One  analytical  development  by  Arnush  [8] 
utilized  radiative  transfer  theory  and  the  small  angle 
approximation  to  characterize  the  light.   Arnush  made  two 
assumptions  that  greatly  weakened  his  model:   (1)  He  assumed 
the  incident  signal  did  not  experience  pulse  broadening  while 
traversing   an  optically  thick  medium,  and  (2)  he  assumed 
the  incident  beam  never  directly  created  internal  emission 
sources.   Stotts  [10]  does  consider  the  previous  details  but 
still  makes  the  small  angle  approximation,  claiming  that  the 
phase  function  is  highly  peaked  in  the  forward  direction. 
Other  attempts  at  gathering  knowledge  include  that  by  Ishimura 
and  Hong  [12]  who  reported  an  analytical  study  of  coherence 
using  first  order  approximations.   The  results  are  valid  for 
weak  fluctuations  in  the  medium. 

Zachor  [16]  uses  a  double-integral  transform  method  which 
is  evaluated  recursively  to  obtain  the  aurole  radiances  contri- 
buted by  successive  scattering  orders.   He  has  assumed  in  his 
calculations  a  homogeneous  unbounded  atmosphere  and  his  results 
are  good  for  short  ranges  only. 

In  any  case,  because  of  the  complexity  of  the  integrals 
involved  ,  exact  analytical  methods  yield  results  only  after 
repeated  numerical  integration.   Even  in  the  small  angle 
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approximation,  the  closed  form  solution  [27]  requires  seven 
successive  integrations  to  obtain  a  single  value  of  the 
irradiance  caused  by  a  unidirectional  point  source  of  light. 
All  this  numerical  work,  besides  being  tedious,  tends  to  mask 
the  functional  dependence  of  the  results  on  the  underlying 
physical  geometry  and  optical  parameters. 

There  are  many  complications  in  solving  the  multiple 
scattering  problem  analytically.   Nevertheless,  analytical 
work  does  contribute  significantly  to  overall  knowledge  of 
the  problem  because  in  many  cases  the  approximations  made  do 
not  substantially  deviate  from  reality. 

In  this  thesis,  two  of  the  more  straightforward  analyti- 
cal methods  are  used;  Gordon  [9]  and  Stotts  [26].  The  former 
concerns  spatial  spread  and  the  latter,  time  spread. 

Dell-Imagine  [15]  derives  a  solution  to  the  radiative 
transfer  equation  analytically  in  the  usual  fashion  using  no 
assumptions  until  he  has  to  actually  get  numerical  results. 
His  numerical  approach  to  evaluating  the  transfer  equation 
seems  quite  attractive  and  is  briefly  mentioned  in  the  follow- 
ing section. 

F.   NUMERICAL  METHODS 

The  solution  of  the  radiative  transfer  equation  is  suffi- 
cient to  specify  the  properties  of  a  received  signal  which  has 
passed  through  a  multiple  scattering  region.   The  solution, 
however,  requires  numerical  computation   on  a  digital  computer. 
The  solution  has  the  general  form  [15], 
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I(x,y,z,t,0 ,0)  =  I(r ,9 ,0,O)exp 

,  t^TT^  2TT/#      ,.        ^-j- 


-Qs 


^JoJoJonir) 


exp 


-Q  c/"n(r)dY 


/  n(r)dX 
!(?**  ,8'  ,0',A) 


P(9,0;  0 ' ,0' )sin0'd0 'd0'dA 


(11) 


where , 


r  =  {(x-nxct) , (y-nyct) ,(z-nzct)} 


(12) 


r'  =  {(x-axc(t-y)) , (y-nyc(t-Y)) , (z-fi  cCt-Y ) ) } 

r"=  {(x-fl  c(t-X)),(y-S2  c(t-X)),(z-Q  c(t-X))} 

fi   =  sin9cos0    Q,      =  sin8sin0    0   =  cos8 
x  y  z 

Q   =  scattering  efficiency  factor    c  =  speed  of  light 

which  has  not  been  restricted  to  any  shape  of  cloud.   By 
establishing  boundary  conditions,  initial  conditions  and 
density  of  particles  n,  a  solution  for  I(x ,y , z , t ,8  ,0)  can  be 
obtained  by  approximating  the  equation  by  a  group  of  simul- 
taneous algebraic  equations.   Dell-Imagine  describes  this 
approximation  and  the  numerical  methods  used  in  solving  the 
equation  and  the  reader  is  referred  to  his  work  for  further 
details . 

G.   MONTE  CARLO  METHODS 

The  Monte  Carlo  method  can  be  applied  to  any  problem  if 
one  knows  the  probability  for  each  step  in  the  sequence  of 


31 


events  and  desires  the  probability  of  the  total  of  all  possi- 
ble events.   Thus,  the  Monte  Carlo  method  may  be  used  to 
study  problems  in  radiative  transfer.   As  mentioned  in  the 
introduction,  many  authors  have  described  their  attempts  to 
do  so.   The  Monte  Carlo  calculations  are  relatively  easier  to 
model  than  other  methods  especially  when  the  geometry  becomes 
difficult.   The  Monte  Carlo  method  is  very  consumptive  of 
computer  time  and  only  approximate  information  is  obtained. 
Some  generally  useful  results  have  been  presented  by  Bucher 
[1] ,  Junge  [6]  and  Hansen  [11].   Monte  Carlo  results  generally 
yield  excellent  agreement  with  experimental  data",  however, 
the  results  deteriorate  statistically  at  large  ranges  or 
narrow  observation  angles. 

Details  of  a  program  for  generation  and  curve  fitting  of 
Monte  Carlo  data  are  contained  in  Junge  [6].   This  thesis 
uses  the  groundwork  of  that  report  to  extend  study  into 
regions  other  than  ultraviolet  light  in  a  homogeneous  space. 

H.   THEORETICAL  RESULTS 

The  purpose  of  this  section  is  to  summarize  the  effects 
scattering  has  on  light  propagation  through  clouds  as  pre- 
dicted by  the  theories  of  the  previous  sections.   In  doing 
this  the  definitions  of  effects  investigated  are  given. 
Actual  quantitative  relations  used  in  this  thesis  are  reserved 
for  later  sections  and  qualitative  descriptions  are  found 
here . 
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Previous  Monte  Carlo,  analytic,  and  numeric  studies  of 
the  multiple  scattering  problem  have  found  that  light  pulses 
are  distorted  in  the  following  manner: 

1.  Angular  Spreading 

Angular  spreading  constitutes  an  effective  decollima- 
tion  of  the  incident  radiation.   It  contributes  to  beam 
spread,  dispersion  in  angle-of-arrival  and  loss  of  spatial 
coherence . 

2 .  Spatial  Spreading 

Spatial  spreading  indicates  the  dimensional  increase 
of  the  beam's  finite  cross  section.   It  is  related  to  the 
above  parameter  in  that  angular  spread  inside  the  medium 
produces  spatial  spread  as  the  pulse  propagates  on.   Spatial 
coherence  and  beam  spread  are  related  to  this  parameter. 

3 .  Multipath  Time  Spreading 

Different  distances  along  various  possible  propagation 
paths  imply  different   transit   times  for  a  photon.   Thus,  a 
short  optical  pulse  will  incur  pulse  broadening  after  tra- 
versing a  multiple  scattering  region.   This  pulse  broadening 
is  called  multipath  time  spreading.   The  maximum  pulse  frequ- 
ency that  can  be  used  in  communicating  is  determined  by  this 
parameter . 

Multipath  time  spread  is  often  defined  as  the  differ- 
ence between  the  average   transit   time  incurred  from  multiple 
scattering  and  the  normal   transit   time  in  the  absence  of 
scattering.   This  definition  is  normally  used  when  the  time 
dependence  of  the  input  is  that  of  a  delta  function.   The 
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previous  description  of  multipath  time  spread  concerns  input 
pulses  of  finite  width. 

4 .   Total  Transmission 

Total  transmission  describes  the  amount  of  irradiance 
left  in  the  pulse  after   traversing   the  medium.   It  is  in- 
versely  related  to  the  beam  attenuation. 

Bucher  [1]  found  that  the  amount  and  distribution 
of  multipath  time  spreading  was  essentially  independent  of 
the  detailed  shape  of  the  scattering  function  for  sufficiently 
thick  clouds.   He  also  observed  that  the  propagation  para- 
meters for  sufficiently  thin  clouds  were  dependent  both  on  the 
cloud  parameters  and  on  the  scattering  function. 

Gordon  [9]  found  that  within  certain  ranges  and  angles 
of  practical  interest,  the  flux  and  beam  spread  function  can 
be   adequately   approximated  by  closed  form  expressions.   His 
work  was  related  to  scattering  under  water  but  can  be  applied 
equally  well  to  atmospheric  scattering  when  the  phase  function 
is  very  forward  peaked. 

Time  spreads  on  the  order  of  microseconds  to  milli- 
seconds have  been  reported  by  Bucher  [1],  Stotts  [10,26]  and 
Ishimura  and  Hong  [12].   Danielson,  Moore  and  van  de  Hulst  [7] 
and  Hansen  [11]  found  that  the  Henyey-Greenstein  function 
adequately  approximated  the  true  cloud  or  haze  phase  function 
in  determining  the  reflection  and  transmission  characteristics 
of  clouds . 

It  is  the  purpose  of  this  thesis  to  further  verify 
results  found  previously  and  to  quantify  some  of  the  more 
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important  cloud  characteristics  with  respect  to  time  and 
spatial  spread. 
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III.   STATEMENT  OF  THE  PROBLEMS 

A.  INTRODUCTION 

There  are  many  areas  in  the  study  of  multiply  scattered 
light  that  could  easily  be  subject  to  further  study.   It  was 
necessary  to  select  a  few  problems  which  were  assailable  using 
available  groundwork  and  techniques. 

B.  THE  PROBLEMS 

The  Monte  Carlo  routine  of  Ref.  6  was  employed  to  investi- 
gate the  following  problems : 

1.  Dependence  of  Spatial  Spread  on  Phase  Function 
Characterize  the  spatial  spread  of  light  traversing 

scattering  media  of  various  types.   Select  phase  functions  of 
high,  moderate  and  low  forward  peakedness  all  with  the  same 
s(cos8>  and  the  same  root  mean  squared  scatter  angle.   Deter- 
mine the  effect  on  spatial  character  of  light  of   using  the 
different  phase  functions.   Determine  the  conditions  concern- 
ing dependence  on  phase  function. 

2 .  Dependence  of  Time  Spread  on  Phase  Function 
Characterize  the  time  spread  of  light  traversing 

scattering  media  of  various  types.   Use  the  same  phase  functions 
selected  above  to  determine  the  effect  on  temporal  character  of 
light   of   using  the  different  phase  functions.   Determine  the 
conditions  concerning  dependence  on  phase  function. 
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3 .  Transition  Region  from  Forward  to  Multiple  Scattering 
Using  the  results  of  the  previous  two  characterizations , 

determine  information  concerning  the  transition  from  the  region 
of  primarily  forward  scatter  to  the  diffusion  region  where 
scatter  is  directed  in  all  directions.   Estimate  the  accuracy 
of  the  information  and  verify  its  agreement  with  other  theories. 

4.  Spatial  Characterization  of  Light  Transversing  Finite 
Clouds 

Characterize  the  spatial  spread  of  light  traversing 

a  homogeneous  medium.   Characterize  the  spatial  spread  of  light 

traversing    the  same  medium  except  that  a  dense  cloud  of 

scatterers  of  finite  thickness  has  been  added.   Compare  the 

two  characterizations  and  discuss  the  results. 

5 .  Model  the  Above  Problems 

Create  a  model  to  simulate  the  geometries  of  the  above 
problems  and  verify  at  each  step  that  the  model  is  indeed 
generating  results  consistent  with  existing  theory. 
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IV.   METHODS  OF  SIMULATION 

A.   MONTE  CARLO  METHOD 

The  computer  routine  of  Ref.  6  was  available  for  use  at 
the  outset  of  this  endeavor.   The  routine  could  calculate 
both  spatial  and  temporal  characteristics  of  light  travers- 
ing a  scattering  medium.   It  required  as  input  the  scattering 
parameters  and  phase  function  parameters  of  the  scattering 
particles  of  the  medium.   The  routine  was  restricted  to  homo- 
geneous atmospheres  utilizing  a  Henyey-Greenstein  phase 
function  or  a  Modified  Henyey-Greenstein  phase  function  [16], 
and  a  given  portion  of  Rayleigh  scattering.   It  was  necessary 
to  modify  the  routine  so  that  it  adequately  represented  the 
geometry  of  a  finite  cloud  and  so  that  any  arbitrary  phase 
function  could  be  adapted  to  it.   These  changes  are  summarized 
in  the  following  sub-sections. 

1.   The  Model 

Figure  4  represents  the  model  used  in  this  simulation. 
The  model  is  identical  to  that  of  Ref.  6  with  respect  to  photon 
path  generation  and  accountability.   The  figure  depicts  a 
cloud  on  whose  left  boundary  the  photons  are  incident  at  the 
center  of  the  accountability  shells,  and  at  whose  right 
boundary  time  and  spatial  information  is  tabulated.   Figure  5 
expands  the  cloud  and  lists  the  names  of  the  parameters  used 
in  the  simulation.   KSCA1  and  KSCA2  represent  the  scattering 
coefficients  in  inverse  kilometers  of  the  inside  and  outside 
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Monte  Carlo  Model 
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Figure  5 
Cloud  Model 


40 


media  respectively.   Likewise,  KABS1  and  KABS2  represent  the 
absorption  coefficients  of  the  two  media.   RPT1  and  RPT2 
represent  the  ratio  of  particulate  to  total  scatter  and  Gl 
and  G2  represent  the  Henyey-Greenstein  phase  function  para- 
meters when  used  in  this  simulation.   THICK  is  the  physical 
thickness  of  the  cloud  in  kilometers.   Appendix  C  gives  details 
on  the  actual  modification  of  the  routine. 

At  each  crossing  of  a  photon  from  inside  of  the  cloud 
to  beyond  the  cloud,  the  total  distance  traveled  is  determined 
and  tabulated  in  time  of  arrival  bins.   Time  spread  and  delay 
information  is  presented  using  these  bins.   A  running  summa- 
tion and  counter  are  used  to  calculate  the  mean  and  standard 
deviation  of  time  required  to  reach  the  cloud  boundary. 

One  other  change  from  the  original  routine  is  conversion 
of  dimensions  from  units  of  extinction  lengths  to  units  of 
kilometers .   Only  minor  changes  in  program  logic  were  required 
to  make  this  change . 

Typical  values  of  atmospheric  parameters  were  used  in 
most  simulations.   A  summary  of  the  parameters  used  is  given 
in  Table  I.   All  diagrams  of  spatial  and  temporal  character 
include  parameters  used  in  the  specific  case. 

The  Monte  Carlo  results  of  Bucher  [1]  were  very  useful 
as  guidelines  for  determining  correct  operation  of  the  model 
created.   Comparisons  of  the  two  models  appear  later  in  this 
section. 
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TABLE  I 
Scattering  Coefficients  used  in  the  Simulation 

KSCACkm"1)  Albedo 

Clear  Atmosphere  .391  1.0 

Light  Haze  .90  1.0 

NOSC  Fog  6.37  1.0 

Water  Cloud  C.2  11.18  1.0 

2 .   Phase  Functions 

As  mentioned  in  the  statement  of  the  problem,  three 
phase  functions  of  different  peakedness  yet  equal  in  ^cos9^ 
were  needed  to  investigate  the  problem.   The  following  phase 
functions  were  employed: 

a.  Henyey-Greenstein  Phase  Function 

Figure  1  depicts  the  functional  dependence  of 
scatter  on  angle  off  direction  of  incidence,  9.   The  closed 
form  expression  for  the  Henyey-Greenstein  function  has  been 
stated  earlier  in  Equation  8.   AG  value  of  .83  was  selected 
so  as  to  match  other  phase  functions  used  in  the  simulation. 
The  peak  of  the  phase  function  is  only  five  units.   Reference 
6  explains  its  use  in  the  Monte  Carlo  routine. 

b.  Numerical  Data  Input  for  Arbitrary  Phase  Functions 
The  other  two  phase  functions  required  development 

of  a  method  for  inputting  phase  functions  consisting  only  of 
data  pairs  and  no  closed  functional  form.   They  are  that 
representing  Water  Cloud  C.2  [21]  at  .53  microns  and  that  of 
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NOSC  Fog.   The  first  was  generated  by  use  of  the  program  of 
Appendix  B,  the  second  by  a  similar  calculation  using  an 
experimentally  obtained  water  particle  distribution.   Appendix 
A  describes  the  details  of  generating  a  random  scatter  angle 
weighted  by  these  arbitrary  phase  functions  defined  by  data 
pairs . 

c.   Rayleigh  Phase  Function 

This  routine  can  also  simulate  a  fractional  amount 
of  the  total  scatter  as  Rayleigh  scatter.   In  this  work  the  use 
of  this  capacity  was  minimal  because  in  most  cases  water  clouds 
of  predominately  large  particles  were  simulated. 
3 .   Automation  of  Contour  Plotting 

Contours  of  equal  relative  photon  flux,  as  described 
in  Ref.  6,  were  used  in  this  work  for  spatial  characterization. 
A  method  of  machine  computation  was  adapted  to  the  DRLITE 
routine.   This  allowed  the  computer  to  perform  the  tedious 
interpolation  made  before  by  pocket  calculator.   Although 
useful  in  some  instances,  the  procedure  proved  quite  computer 
time  consumptive.   Its  usefulness  was  very  pronounced  in  comparing 
the  Monte  Carlo  routine  with  Gordon  [9]  as  will  be  discussed 
later.   Appendix  D  gives  details  of  the  adaption  of  this  method 
to  the  routine. 

B.   ANALYTICAL  METHODS 

As  explained  in  the  previous  section  on  theoretical  pre- 
liminaries, analytical  methods  ordinarily  lead  to  solutions 
which  are  applicable  only  under  certain  conditions.   This  is  so 
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because  in  order  to  reach  a  usable  solution,  approximations 
have  to  be  made.   Nevertheless,  in  the  regions  of  their 
applicability,  analytical  solutions  are  excellent  sources  of 
the  information  needed  to  verify  results  generated  by  other 
more  universally  applicable  methods .   In  this  work  two  analyti- 
cal treatments  were  employed  for  just  that  reason.   Gordon's 
[9]  paper  on  practical  approaches  to  underwater  scattering 
was  used  to  verify  the  spatial  output  of  the  Monte  Carlo  routine 
created  and  Stott's  [26]  closed  form  expression  for  time 
spread  was  used  to  verify  the  temporal  output  of  the  program. 
The  two  treatments  are  summarized  in  the  following  paragraphs. 
1.   Effective  Attenuation  Coefficient  (EAC)  Method 

Gordon's  paper  [9]  presents  the  concept  of  an  effective 
attenuation  coefficient  which  considerably  simplifies  under- 
water multiple  scattering  calculations.   Concise  expressions 
for  the  total  flux  through  an  aperture  and  the  beam-spread 
function  are  defined  in  terms  of  the  EAC.   The  derivation  of 
these  expressions  is  included  for  easy  reference  as  Appendix  E 
as  is  a  documentation  and  description  of  the  program  adapting 
them  to  machine  computation.   It  is  interesting  to  note  that 
Gordon's  treatment  is  actually  a  truncated  version  of  the 
exact  analytical  solution  of  the  multiple  scattering  problem 
using  the  small  angle  approximation  provided  by  him  [28].   The 
truncation,  in  effect,  considers  the  scattering  phase  function 
as  a  forward  scattering  delta  function,  thus  does  not  consider 
higher  order  terms  of  the  exact  solution.   It  is  not  surprising 
then  that  the  result  is  only  applicable  under  certain  conditions. 
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Fortunately  these  conditions  are  found  often  in  the  underwater 
optical  world.   The  comparison  of  this  method  to  the  Monte 
Carlo  method  is  found  in  the  next  section.   It  should  be  kept 
in  mind  that  Gordon's  treatment,  even  in  the  truncated  manner, 
closely  agrees  with  experimental  data  [29].   The  reader  is 
referred  to  Appendix  E  for  further  details  on  the  EAC  method. 
2 .   Closed  Form  Time  Spread  Expression 

Stott's  treatment  [26]  of  time  spread  in  a  multiple 
scattering  region  is  short,  concise  and  seems  to  agree  well 
with  existing  time  spread  data.   It  was  therefore  chosen  as 
an  additional  means  for  verification  of  time  spread  results 
given  by  the  Monte  Carlo  method.   The  details  in  deriving 
Stott's  closed  form  solution  for  time  spread  are  located  in 
Appendix  G  for  easy  reference.   The  result  is  the  average 
additional  multipath  time  required  for  a  photon  to  traverse 
a  scattering  medium  with  physical  thickness  z,  optical  thick- 
ness t,  albedo  aQ   and  RMS  scatter  angle  yQ.      The  time  spread 
is 


L  =  z  (_^  C(l+f«0TY0f)1- •-!]-!)      (1 

C  \C0oTYoZ  / 


2) 


This  result  is  compared  to  Monte  Carlo  results  in  later  sec- 
tions . 

C.   COMPARISON  OF  METHODS  EMPLOYED 

A  problem  usually  encountered  in  mathematically  modeling 
any  physical  situation  is  to  ensure  at  all  times  that  the  model 
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employed  is  turning  out  reasonable  answers ,   To  verify  that 
this  was  in  fact  the  case,  frequent  cross-checks  were  made 
between  the  models.   The  proof  of  the  theory  follows,  of 
course,  in  agreement  between  the  models  and/or  agreement  with 
experimental  observations.   This  section  compares  the  results 
of  the  Monte  Carlo  model  with  the  analytical  models  as  well 
as  with  the  results  of  Bucher ' s  [1]  independently  created 
Monte  Carlo  model.   Once  confidence  is  established,  conclu- 
sions can  be  drawn  about  the  regions  of  effectiveness  of  each 
model.   Other  more  routine  checks  for  error  are  noted  in 
Appendix  F. 

1.   Monte  Carlo  vs.  Analytical 
a.   Spatial  Characterization 

The  EAC  method  derived  and  verified  in  Appendix  E 
is  compared  to  the  Monte  Carlo  method  in  this  section.   First, 
the  method  of  spatial  characterization  adapted  from  Ref.  6  is 
briefly  explained. 

(1)  Description  of  Relative  Flux  Contours .   At 
each  angular  bin  of  each  reference  shell  the  relative  flux  is 
given  by 

_  ,  .  .    „,         #  Photons  Passing  Through  Bin     , ,  „ , 

Relative  Flux  =  =  ,  -, — Tj — r- =n qr— § a 5  5  !    (13) 

Total  Number  Ejected  x  Area  of  Bin 

The  negative  logarithm  of  this  quantity  is  calculated  for  each 
bin.   This  array  of  values  is  interpolated  as  described  in 
Appendix  D, and  contours  of  equal  relative  photon  flux  are 
drawn.   Likewise,  the  negative  logarithm  of  the  beam  spread 
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function  evaluated  by  the  EAC  method  is  evaluated  at  a  similar 
array  of  spatial  positions  and  equal  flux  contours  are  drawn. 
The  two  sets  of  contour  lines  are  compared  in  the  following 
paragraphs . 

(2)  Comparison.   Flux  contours  calculated  using 
Henyey-Greenstein  G  parameters  of  .70  and  .95  for  each  albedo 
of  .80  and  .90  were  drawn  at  integer  values  of  the  negative 
logarithm.   Figures  6  through  9  show  these  comparisons.   The 
dotted  flux  contours  are  those  of  the  EAC  method  from  zero  to 
45  degrees  and  the  solid  contours  are  those  of  the  Monte  Carlo 
method.   The  term  S  on  the  graph  is  the  scattering  coefficient 
as  well  as  the  albedo  since  the  extinction  coefficient  is 
unity  in  each  case.   Dimensions  are  in  kilometers  but ■ in  this 
case  one  kilometer  is  equivalent  to  one  extinction  length. 

In  general,  agreement  is  good  between  the 
methods  for  highly  peaked,  G  =  .95,  phase  functions  at  dis- 
tances less  than  10  extinction  lengths  and  angles  varying  from 
one  to  20  degrees.   This  comparison  agrees  with  results  given 
in  Gordon  [9].   It  should  be  noted  that  angular  agreement  is 
nearly  independent  of  albedo  whereas  range  agreement  is  better 
when  large  absorption  is  present. 

(3)  Conclusions .   Inside  of  a  few  degrees  at  moder- 
ate to  large  extinction  lengths  the  EAC  method  becomes 
undependable  in  most  cases .   The  Monte  Carlo  and  EAC  methods 
agree  well  in  regions  where  agreement  should  be  expected.   The 
Monte  Carlo  method  has  the  advantage  over  the  EAC  method  when 
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considering  atmospheric  scattering  because  in  this  case 
absorption  is  small  and  large  optical  thicknesses  are  encoun- 
tered often.   Confidence  can  be  placed  in  the  Monte  Carlo 
routine  for  spatial  characterization  of  light  in  a  multiple 
scattering  region. 

b.   Temporal  Characterization 

Stotts's  [26]  closed  form  expression  for  time  spread 
which  is  summarized  in  Appendix  G  is  compared  to  time  spread 
found  using  the  Monte  Carlo  method  in  this  section.   The 
technique  used  in  tabulating  time  spread  information  in  the 
Monte  Carlo  model  was  explained  earlier  in  the  section  titled 
"The  Model". 

(1)  Comparison .   Multipath  time  spread  is  defined 
as  the  difference  between  the  average   transit   time  incurred 
from  multiple  scattering  and  the  normal   transit   time  in  the 
absence  of  scattering.   The  expression  for  time  spread  derived 
by  Stotts   is  given  in  Equation  12.   This  equation  is  plotted 
in  Figures  10  through  14  with  an  albedo  of  one,  an  RMS  scatter 
angle  of  .65  and  optical  thickness,  x,  ranging  from  zero  to 
100  for  physical  thickness  values,  z,  of  .5,  1.0,  2.0,  3.0  and 
1.0  kilometers.   Monte  Carlo  data  points  are  marked  with  X's 
using  various t  values  with  other  parameters  the  same.   Also 
plotted  are  Bucher's  [l]  data  which  is  explained  later.   In 
general,  agreement  is  within  one  microsecond  for  optical  thick- 
nesses less  than  10  and  within  five  microseconds  for  optical 
thicknesses  between  10  and  60. 
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(2)  Conclusions.   The  Monte  Carlo  and  closed  fo 


rm 


methods  agree  well  in  regions  where  agreement  should  be  expec- 
;j  ted.   Because  of  this  agreement,  confidence  can  be  placed  in 

the  Monte  Carlo  routine  for  temporal  characterization.   The 
, Monte  Carlo  routine  has  the  advantage  over  the  closed  form 
"expression  in  that  different  phase  functions  can  be  simulated, 

thus  allowing  study  of  the  effects  of  the  details  of  the  phase 

functions  on  time  spread  of  light  in  a  scattering  medium. 

Because  time  spread  is  expressed  quite  well  by  the  closed  form 

expression,  little  dependence  on  these  details  is  expected. 

The  dependence  is  discussed  in  later  sections. 

2.   This  Monte  Carlo  Model  vs.  Bucher's  Monte  Carlo  Model  [l] 


a.   Spatial  Characterization 

The  Monte  Carlo  method  is  tested  for  agreement 
with  Bucher's  [1]  Monte  Carlo  results.   Because  both  methods 
are  built  with  nearly  the  same  geometry,  nearly  identical 
results  are  expected.   First  the  graphic  representation  of 
data  used  by  Bucher  needs  explanation. 

(1)  Description  of  Scaling  Method.   Bucher  [1] 
defines  for  a  scattering  medium  the  diffusion  distance 


D.  =  t— J2 — _-   ,  (14) 

d   1-<cos0> 

where  D  is  the  mean  free  path  between  collisions  and  <cos9)>  is 
the  average  cosine  of  the  scattering  angle  as  before.   Using 
this  formula,  the  effective  scattering  thickness  t  ,  is  given  by 
the  relation 


58 


Td  =  D^  =  T(l-<cose>)  (15) 

d 

where  T  is  the  physical  thickness  of  the  cloud  and  x    is  the 
usual  optical  thickness.   The  spatial  results  are  then  pre- 
sented by  graphically  displaying  the  average  displacement 
<r>  in  relation  to  the  diffusion  thickness,  D , ,  as  a  function 
of  effective  scattering  thickness  t,.   Also  presented  is  the 
I  critical  transverse  displacement,  r  ,  in  relation  to  D,  inside 
which  50  per  cent  of  the  photons  exit  the  cloud  at  effective 
scattering  thickness,  x,.   The  reader  is  referred  to  Bucher 
[1]  for  further  details. 

(2)  Comparison.   Monte  Carlo  calculations  were 
made  and  interpreted  in  terms  of  the  units  described  above. 
Figures  15  and  16  show  the  results  at  various  effective 
scattering  thicknesses  using  an  albedo  of  unity  and  a  Henyey- 
Greenstein  function  for  which  G  =  .83.   The  X's  represent  data 
collected  here  and  the  continuous  line  is  a  best  fit  curve  to 
Bucher 's  data  using  the  same  parameters.   As  can  be  seen, 
agreement  is  excellent  in  both  cases. 

(3)  Conclusion.  Spatial  characterization  by  this  Monte 
Carlo  routine  is  consistent  with  a  previously  developed  routine. 
This  creates  confidence  in  the  routine.   It  can  now  be  used  to 
study  the  effect  phase  function  details  have  on  the  spatial  char- 
acter of  light. 

b.   Temporal  Characterization 

The  time  soread  data  of  the  Monte  Carlo  routine  is 
compared  to  Bucher  [1]  in  this  section.   Two  methods  of 
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presentation  are  used.   One  uses  time  bins,  the  other  uses 
Bucher's  best  fit  equation  for  time  spread. 

(1)  Comparison.   Bucher's  equation  for  multipath 
time  spread  as  defined  in  an  earlier  section  is, 

e  9  rn        94  10  94 

L(tJ  =  i2£±(T  )•   =  3.9i  x  io_   T(t  J"   (16) 

d     c    d  d 

where  L(t,)  is  the  time  spread  in  seconds,  T  is  the  physical 
thickness  of  the  cloud  in  meters,  x,  is  the  effective  scatter- 
ing thickness  of  the  cloud  and  c  is  the  speed  of  light  in 
meters  per  second.   This  function  is  plotted  alongside  Monte 
Carlo  data  and  Stott ' s  closed  form  equation  in  Figures  10 
through  14  for  an  albedo  of  unity,  a  Henyey-Greenstein  phase 
function  with  G  =  .83,  and  a  physical  thickness  of  .5,  1.0, 
2.0,  3.0  and  4.0  kilometers. 

Agreement  is  well  within  one  microsecond  for 
nearly  the  entire  range  which  is  well  within  statistical 
error.   In  these  trials  only  50  -  500  photon  histories  were 
tabulated  which  seems  to  be  an  indication  that  collection  of 
time  spread  information  does  not  require  unwieldy  amounts  of 
computer  time . 

Another  comparison  of  time  spread  information 
is  graphically  displayed  in  Figures  17  and  18.   The  number  of 
photons  arriving  in  each  time  bin  is  normalized  to  the  maximum 
number  in  any  bin  and  plotted  as  a  function  of  time  in  micro- 
seconds.  Bucher's  data  is  superimposed  by  a  dotted  line  on 
the  time  bin  data  presented.   Figure  17  plots  results  for  an 
optical  thickness  of  15  and  Figure  18  for  an  optical  thickness 
of  30.   Other  parameters  are  listed  on  the  graphs.   Agreement 
once  again  is  very  good  between  the  methods. 
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(2)  Conclusion.   Once  again,  confidence  is  ensured 
in  temporal  characterization  of  time  spread  and  pulse  spread 
by  the  Monte  Carlo  routine.   The  routine  can  now  be  used  in 
investigating  time  spread  dependence  on  phase  funtion. 
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V.   RESULTS 

A.   SPATIAL  SPREAD  DEPENDENCE  ON  PHASE  FUNCTION 

In  this  section  five  graphs  are  used  to  display  the 
dependence  of  the  spatial  spread  of  light  on  phase  function 
when   traversing   a  scattering  medium.   All  five  diagrams 
show  the  relative  flux  contours  described  earlier  in  this 
thesis  and  in  Ref.  6.   Each  figure  lists  the  parameters  and 
the  phase  function  used  in  the  simulation.   Each  axis  is  in 
units  of  kilometers.   Incidence  of  photons  is  at  the  origin 
and  to  the  right  in  all  cases.   The  left  edge  of  the  cloud, 
if  it  exists ,  is  along  the  axis  perpendicular  to  direction 
of  incidence  and  the  right  edge  of  the  cloud  is  indicated  by 
a  dotted  line  if  it  is  within  the  boundaries  of  the  graph. 
For  easy  reference,  Table  II   tabulates  the  parameters  and 
phase  function  used  in  Figures  19  through  36. 

The  phase  function  of  Water  Cloud  C.2  is  used  in  Figure  19 
to  show  the  peakedness  of  the  relative  flux  contours  out  to 
at  least  20  extinction  lengths  when  the  albedo  is  only  .90. 
A  "toe"  appears  on  the  contours  at  small  angles  due  to  the 
high  degree  of  forward  scattering  of  Water  Cloud  C.2. 

Figures  2  0  and  21  compare  the  flux  contours  due  to  NOSC 
Fog  and  Henyey-Greenstein  using  G  =  .83,  respectively.   The 
flux  contour  for  NOSC  Fog  is  more  forward  peaked  than  that  of 
the  Henyey-Greenstein  function  for  distances  less  than  three 
kilometers  or  19  extinction  lengths  but  flux  contours  are  very 
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similar  at  distances  greater  than  three  kilometers.   Outside 
of  2  0  extinction  lengths  the  contours  and  therefore  the 
spatial  character  of  the  light  they  represent  become  indepen- 
dent of  the  phase  function  used.   The  contours  of  these  two 
graphs  are  much  less  forward  peaked  than  that  of  Figure  19 
because  the  albedo  in  the  former  is  unity. 

Figures  2  2  and  2  3  show  similar  results  when  comparing 
results  when  Water  Cloud  C.2  [21]  and  Henyey-Greenstein  were 
used,  respectively.   In  all  figures  the  cloud  was  extended  to 
infinity  at  the  right  of  the  origin,  explaining  the  discontinu- 
ity in  contour  slope. 

B.   TEMPORAL  SPREAD  DEPENDENCE  ON  PHASE  FUNCTION 

In  this  section  nine  histograms  and  one  table  are  used  to 
display  the  dependence  of  time  spread  of  light  on  phase  func- 
tion when   traversing  a  scattering  medium.   All  nine  diagrams 
list  the  parameters  in  the  same  fashion  as  in  the  last  section. 
The  essential  parameters  can  also  be  found  in  tabulated  form 
in  Table  II.   The  information  is  plotted  as  the  number  of 
photons  arriving  in  any  bin  normalized  to  the  bin  of  most 
frequent  arrival  versus  the  time  spread  incurred.   Note  that 
the  bins  are  logarithmic  and  therefore  represent  much  less 
span  in  time  to  the  left  than  to  the  right. 

Figures  24-,  25  and  26  compare  time  spread  for  NOSC  Fog, 
Water  Cloud  C.2  and  Henyey-Greenstein,  respectively  in  a  cloud 
one  kilometer  thick  of  optical  thickness  four.   It  can  be  seen 
by  comparing  the  graphs  that  at  this  optical  thickness  the 
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higher  peaked  phase  functions  have  more  photons  arriving  at 
earlier  times. 

Figures  27,  2  8  and  2  9  similarly  compare  time  data  for  a 
!  cloud  with  a  physical  thickness  of  one  kilometer  and  an  optical 
1  thickness  of  eight.   The  time  spread  diagrams  begin  to  look 
very  much  alike  but  still  vary  slightly  with  peakedness.   At 
this  thickness  the  average   transit   time  is  nearly  identical 
for  the  three  cases  but  there  are  still  photons  of  the  higher 
peaked  phase  functions  which   traverse   the  cloud  with  very 
little  time  delay.   This  indicates  some  dependence  on  phase 
function  for  this  thickness. 

Figures  30,  31  and  32  conclude  the  graphical  time  spread 
presentation  with  a  cloud  one  kilometer  thick  with  an  optical 
thickness  of  15.   There  is  no  noticeable  difference  in  the 
graphs  that  can  not  be  attributed  to  the  statistical  deviation 
expected  in  a  Monte  Carlo  routine.   Table  III  summarizes  the 
time  spread  data. 

From  these  results  it  can  be  said  that  the  time  spread 
becomes  independent  of  the  details  of  the  phase  function  at 
optical  thicknesses  greater  than  15. 

C.   REGION  OF  TRANSITION  FROM  FORWARD  TO  MULTIPLE  SCATTER 

The  region  of  transition  from  forward  to  multiple  scatter 
can  be  defined  as  that  region  for  which  (1)  distances  of  lesser 
magnitude  show  time  and  spatial  spread  depending  markedly  on  the 
phase  function  of  the  scattering  particles  and  (2)  distances  of 
greater  magnitude  show  that  time  and  spatial  spread  do  not  depend 
on  the  phase  function  of  the  scattering  particles. 
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TABLE  III 

Summary  of  Time  Spreading 
by  One  Kilometer  Thick  Cloud 


KSCA1 

H.G.  G=.83 

W.C.  C.2 

NOSC  Fog 

(km-1) 

(y-sec) 

^  Mean 
Sigma 

1.7 

1.2 

1.2 

3.1 

2.4 

2.8 

0  Mean 
Sigma 

2.9 

2.9 

2.9 

4.2 

4.5 

4.1 

15  Jean 
Sigma 

4.4 

4.6 

4.3 

4.4 

5.  3 

4.7 

30  ^ean 
Sigma 

7.5 

7.1 

7.9 

5.  3 

5.6 

5.5 
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From  the  previous  two  sections  the  region  of  transition 
from  forward  to  multiple  scatter  is  between  about  10  -  20 
optical  thicknesses.   This  agrees  favorably  with  observations 
made  by  other  authors  ["+0]. 

D.   EFFECT  ON  SPATIAL  CHARACTER  OF  LIGHT  WHEN  PASSING  THROUGH 
A  CLOUD 

In  this  section  four  graphs  are  used  to  display  the  effect 
of  finite  thickness  clouds  on  the  spatial  spread  of  light. 
All  four  diagrams  use  the  relative  flux  contours  and  the  list 
of  parameters  defined  in  earlier  sections  of  this  report. 

Figure  3  3  shows  the  relative  flux  contours  of  an  infinite 
very  clear  atmosphere.   Figure  34  shows  the  effect  when  a  cloud 
of  optical  thickness  about  nine  is  placed  in  the  light's  path. 
Deformation  of  the  contours  is  very  obvious. 

Likewise,  Figures  35  and  36  show  the  effect  of  a  . 5  kilo- 
meter Water  Cloud  C.2  when  placed  in  the  path  of  light  in  a 
clear  atmosphere.   The  sharp  peaks  on  the  flux  contours  are 
rounded  and  as  expected,  light  is  significantly  redirected 
when  passing  through  the  cloud.   Clouds  of  large  optical  thick- 
ness would  eventually  cause  the  Lambertian  irradiance  well 
known  in  the  theory  of  light  propagation. 
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VI.   DISCUSSION 

In  this  thesis  Monte  Carlo  methods  and  simple  analytical 
methods  were  used  to  characterize  light  transmission  through 
a  multiple  scattering  medium.   Complicated  numerical  methods 
were  not  used  nor  were  highly  mathematical  models .   This  type 
of  approach  remains  for  future  investigations.   It  is  antici- 
pated that  results  drawn  from  present  methods  will  be  confirmed, 

Much  information  on  the  multiple  scattering  problem  re- 
mains to  be  gathered.   The  present  Monte  Carlo  routine  can  be 
adapted  to  numerous  useful  geometries  and  atmospheres  with 
gradients  and  time  variation  of  physical  parameters .   Propaga- 
tion over  and  around  solid  bodies  such  as  mountains  and  the 
earth  itself  can  be  modeled.   Layered  media  such  as  cloud, 
air  and  water  interfaces  can  also  be  simulated  by  the  routine. 


89 


VII.   CONCLUSIONS 

A  computer  routine  based  on  a  Monte  Carlo  model  was  de- 
veloped to  simulate  light  propagation  through  a  scattering 
absorbing  medium.   The  routine  can  simulate  a  plane  parallel 
cloud  of  scatterers  with  the  desired  parameters  located  within 
a  medium  with  another  set  of  parameters.   The  phase  functions 
of  each  medium  can  be  arbitrarily  defined  using  a  set  of  data 
pairs  or  can  be  approximated  by  closed  form  expressions.   The 
data  generated  by  the  routine  has  been  checked  for  accuracy 
against  other  theories  using  analytical  methods.   Each  compari- 
son shows  adequate  agreement  between  theories  where  agreement 
can  be  expected. 

The  routine  will  automatically  plot  spatial  information 
necessary  in  characterizing  light  transfer  through  the  medium 
by  use  of  contours  of  equal  photon  flux.   It  will  also  generate 
histograms  depicting  time  spread  information  for  light  passing 
through  finite  clouds . 

Using  the  Monte  Carlo  routine  created  and  inputting  phase 
functions  of  different  peakedness,  it  has  been  found  that  both 
spatial  and  temporal  spread  are  independent  of  the  details  of 
the  phase  function  for  thicknesses  greater  than  15  extinction 
lengths.   The  region  of  transition  from  forward  scatter  to 
multiple  scatter  is  between  10  -  20  extinction  lengths. 

The  routine  has  also  been  used  to  study  the  effect  finite 
thickness  clouds  have  on  the  spatial  character  of  light. 
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APPENDIX  A 

I.   GENERATION  OF  A  RANDOM  SCATTER/ ANGLE 
WEIGHTED  BY  AN  ARBITRARY  PHASE  FUNCTION 

A.  STATEMENT  OF  THE  PROBLEM 

Reference  6  outlines  the  theory  and  calculations  necessary 
to  generate  random  numbers  weighted  in  accordance  with  functions 
that  represent  characteristics  of  a  scattering  medium.   At  each 
collision  new  scatter  angles  and  a  distance  were  generated  by 
closed  form  expressions  which  weight  a  uniformly  distributed 
random  number.   However,  in  the  more  general  case  of  a  poly- 
dispersion,  the  computed  phase  function  could  not  be 
represented  adequately  in  closed  form.   This  made  it  impossible 
to  invert  the  function  enabling  analytical  generation  of  a 
weighted  scatter  distribution.   It  was  necessary  to  create  a 
method  for  generation  of  a  random  theta  weighted  by  an  arbitrary 
phase  function. 

B.  METHOD  USED  IN  WEIGHTING  THETA 

Reference  21  and  Appendix  B  explain  in  great  detail  a 
method  for  calculating  the  representative  phase  function  of  a 
polydispersion  given  the  particle  size  distribution  and  com- 
position.  Using  this  computer  adaptation  of  Mie  theory  ,  values 
of  the  averaged  normalized  phase  function  were  generated  at 
selected  angles  of  scatter.   Given  this  data  and  a  random 
number,  R,  weighted  uniformly  over  the  closed  interval  [0,  1], 
the  problem  is  to  solve  the  following  equation  for  9R. 
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f 

r  =  ±   P(9)sin  9  d9  (1) 


I 


7rP(0)sin  9  d0 


'0 

where  P(9)  is  the  averaged  normalized  phase  function  and  9  is 
the  variable  of  integration.   Notice  that, 

9=0  when  R  =  0, 

(2) 
0  =  77  when  R  =  1 , 

as  you  would  expect  on  the  closed  interval.   P(9)  is  not  a 
continuously  defined  function  over  the  interval  [0,  tt  ]  in  this 
case  so  the  integrals  are  represented  numerically  by  the 
trapezoidal  rule.   Figure  37  diagrams  the  basic  procedure  used, 
N  values  of  P(9)  are  selected  so  as  to  closely  approximate  the 
phase  function.   Of  course,  more  values  are  selected  near  the 
small  scatter  angles  to  adequately  describe  the  sharp  peak. 
The  N  values  are  multiplied  by  the  sine  of  the  respective 
scatter  angle  (P(0)  includes  sin  9  implicitly  in  Figure  37). 
N-l  trapezoids  are  established  using  these  N  values.   The  area 
of  each  trapezoid  is  divided  by  the  total  area  of  all  trape- 
zoids thus  establishing  a  weight  for  the  respective  interval. 

(p.sin9 .+p.  ..sine .., )(e:..n-e.) 

T7  IX    1  +  1      1+1     1  +  1    i        f  0  >, 

w.  =vr_i v.3; 

1   v ^1CP.sin9.+P.i1sin9.x1  )(Q. .,-9.) 

2.,  i   i  i+i   i+i   i+i  i 

i=l 
This  of  course  requires  that 


N-l 

S  w    =  i.  ^) 

n=l  " 
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Figure  37 
Diagram  of  Random  Generation  of  Theta 
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With  the  weight  of  each  panel  known,  the  uniform  random 
number  R  is  used  to  solve  for  the  first  value  of  n  such  that, 

n 

£   Wi  >  R  (5) 

i=l 

is  satisfied.   The  random  9-d  required  is  somewhere  between  9 

R   n  n 

and  9  L, .   Within  the  panel,  9  is  weighted  linearly  in  a  fashion 
n+1  e  '  to  J 

as  explained  in  Ref.  6.   For  a  given  panel,  9.,  9.   ,  P.  and 
P.+,  are  known.   From  these  the  slope  and  intercept  are  estab- 
lsihed. 


P.   -P. 

B.  =   X  X   x  A.  =  P.-B.9.       (6) 

l    8  ... -9  .  l     l   x  x 

i+l   l 


The   panel    is    normalized   as    follows : 

/i+l 
NORM/  (A.+B.9)d9    =    1  (7) 


i      l 


so   that 


NORM   =    g-: (8) 


All   that   remains    is    to    solve    the    following    for 

(r-x;  w±) 


R: 


R,    =    — x1"1 =    NORM/    *tA.+B.9)d9         (9) 

1         (nAVr«.) 


£"±-£"1' 

i=l        i=l 


l£\+Bi{ 


9U 


After  a  few  algebraic  steps,  the  solution  is 


^feM# 


+    Bi>0 
-    Bi<0 


(10) 


where , 


:i  ■  Ri{Ai»i*/JTJ)*l1-Rj{^tAih)  (11) 


This  9D  has  the  desired  properties.   The  FORTRAN  coding  of  this 
method  is  found  in  the  Subroutine  RANTH. 
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APPENDIX  B 

I.   DESCRIPTION  AND  DOCUMENTATION  OF  PROGRAM 
TO  ADAPT  MIE  THEORY  TO  MACHINE  COMPUTATION 

A.   INTRODUCTION 

Mie  theory  has  been  adapted  to  machine  computation  on  many 
occasions.   Deirmendjian  [21]  provides  an  excellent  guideline 
for  using  Mie  theory  on  spherical  polydispersions .   Using 
Deirmendjian' s  guideline,  a  computer  routine  was  created  to 
generate  the  volume  scattering  and  extinction  cross  sections 
and  the  corresponding  elements  of  the  normalized  scattering 
matrix  for  a  polydispersion  where  the  number  of  particles  per 
unit  volume,  per  unit  radius  is  given  by, 


n(r)  =  ar  e  o<r<°°  (1) 


where  r  is  the  particle  radius.   The  four  constants  a,  a,  b 
and  y  are  positive  and  real  and  a  is  an  integer.   They  are  not 
independent  of  each  other,  and  are  related  to  quantities  in  the 
particle  frequency  distribution.   The  radius  which  is  most 
frequently  found  in  the  particle  distribution  is  r   and  N  is  the 
total  number  of  particles  per  unit  volume.   Both  N  and  r   can 
be  found  by  experimental  measurement.   In  terms  of  N  and  r   the 
constants  of  the  distribution  are  found  using: 


vf 


-OD         v 

i   a   -br'    _     a     r,a+l.       ,  0  >. 
N  =  a/re      dr  =  K— — )       (2) 

(q  +  1)     Y 

Yb  Y 
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b  =  -J2L- 

Y  (3) 

yr 
'  c 

and  by  choice  of  a.   and  y    to  best  fit  the  experimental  distri- 
bution function.  T    is  the  usual  gamma  function. 

B.   EQUATIONS  TO  BE  CALCULATED 

The  equations  used  to  generate  scattering  data  for  a 
polydispersive  system  are  extensions  of  those  used  in  a  mono- 
dispersion.   The  distribution  has  been  created  in  a  manner 
requiring  that 


A 


N  =  /   n(r)dr  (4) 

where  n(r)  is  a  continuous  and  integrable  function  within  the 
range  and  represents  the  partial  concentration  per  unit  volume 
per  unit  increment  of  radius  r.   The  values  of  interest  are 
volume  scattering  cross  sections  and  the  corresponding  elements 
of  the  normalized  scattering  matrix  [18,  21].   These  values 
can  be  computed  using 

•oo 


/•UU 

TTk"3/       x2n(x] 


3         [A,    n(x)]    =    Trk    3/       x2n(x)Qsca(x)dx         (5) 

S  CcL 


and 


f 

Pext^»    nCx)]    =    TTk    3/       x2n(x)Qext (x)dx         (6) 

^CO  0 

Pj(9)    =    k-IicT/      nU)    ij(0)dx        3*1,2  (7) 


o 
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where 

Q    (x)  =  scattering  efficiency  factor 
sea  °  J 

Q   ^(x)  ~  extinction  efficiency  factor 
i.(x)  =  dimensionless  intensity  parameters 

Each  of  the  above  terms  will  be  examined  in  detail  in  following 
sections . 

C.   CALCULATION  OF  SCATTERING  EFFICIENCY  FACTORS  AND  DIMENSION- 
LESS  INTENSITY  PARAMETERS 

To  perform  numerical  integration,  the  integrand  must  be 
evaluated  at  many  different  values  within  the  range  of  inte- 
gration.  In  this  case  the  scattering  and  extinction  efficiency 
factor  must  be  evaluated  for  numerous  Mie  size  parameters  x. 
Scattering  efficiency  is  a  name  given  to  the  ratio  of  total 
scattering  cross  section  to  geometrical  cross  section.   A 
similar  definition  is  used  to  define  extinction  efficiency. 
For  any  particle  size  and  composition,  Mie  theory  gives  the 
total  scattering  cross  section  as  [18,  21], 


Jn 


6    (m,x)  =  i/   (A,A*+A0A*)d(o  (8) 

sea         2  /     112  2 


where 


00  2n+l 


kA1  =  Sl(m,x,9)  =  £  HTnTIT^nVV^    (9) 


n=l 

GO 


2n+l 


kA2  =  S2(m,x,9)  '-   gi-^i-7(bn7rn+anTn)    (1 


0) 


9  8 


using  the  complex  index  of  refraction  of  the  particle  m,  the 

Mie  coefficients  a  ,  b   and  angular  coefficients  tt  ,  t  ,  each 

n'   n       b  n    n ' 

of  which  is  described  in  detail  in  later  sub-sections.   Using 
these  equations,  the  scattering  efficiency  factors  are  [18,  21], 


9   °° 
Qsca(m,x)  =^E  (2n+l)(|an|2+|bn|2),    (11) 

x  n=l 

Q   .(m,x)  =  =hr   £  (2n+l)Re{a  +b  }.        (12) 
xext  2  ■?-?-,  n   n 

x  n-1 


The  dimensionless  intensity  parameters  are  given  by  the 
expressions , 


i  (x,m,8)  =  k2A1A='t  =  SiS{: 


i2(x,m,9)  =  k2A2A^  =  S2S2* 


(13) 


where  A,,  A? ,  S,  and  S_  are  the  same  as  above.   The  procedure 

used  in  calculating  a  ,  b  ,  t   and  it   is  described  next. 

&   n'   n5   n      n 

1.   Computation  of  t   and  tt 
n n 

The  two  angular  functions  t   and  tt   are  required  before 
6  n      n        H 

the  dimensionless  intensity  parameters  can  be  calculated.   t 

and  it   are  functions  of  u  =  cos0  only  and  as  the  subscript 

implies,  there  are  actually  many  different  angular  functions. 

Both  functions  are  defined  in  terms  of  Legendre  polynomials  [18] 

and  their  derivatives.   The  following  recursion  relations  are 

used  in  calculating  t   and  tt  , 

&   n      n' 
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TT      (0)      S     cOSe^^    TT         _C6) ~    TT         .(6)  (14) 

n  n-1      n-1  n-1      n-2 


t    (9) 
n 


cos 


9  [it    (9  )-tt         (9)]     -  [(2n-l)sin29    tt      .(6)1 
L  n  n-z       JJ       "-  n-1        J 


+    t      9(0)  (15) 

n-z 


with 


tt    (9)    =    0  t    (9)    =    0 

0  o 

tt,  (9)    =    1  t    (9)    =    cos    9  (16) 

1  o 

tt2(9)    =    3    cos9  t„(9)    =    3    cos    20 

2 .   Computation  of  a   and  b 
n n 

The  Mie  coefficients  a   and  b   have  a  complex  depen- 

n       n  t-  t- 

dence  on  the  index  of  refraction  of  the  particle,  and  the 
surrounding  medium  and  also  depend  on  the  size  parameter,  x. 
There  exists  many  different  expressions  for  the  Mie  coeffici- 
ents in  terms  of  known  mathematical  functions.   The  form  used 
in  this  program  is  described  in  terms  of  spherical  Bessel 
functions  each  of  which  is  described  in  later  sections.   The 

relations  for  a   and  b   are 
n      n 


a  = 


b  =  £ 

n   - 


(17) 


u[yjn_1(y)-njn(y)jxjn(x)-eyjn(y)|^xjn_1(x)-njn(x)j 
n  u[yjn^1(y)-njn(y)]xhn2(x)-eyjncy)[hj_1(x)-nh2  (x)J 

|_yjn_1(y)-njn(y)j  xjn(x)-yyjn(y)  |^xjn_1(x)-njn(x)j 
ryJn_1Cy)-nJn(y)]xh^(x)-yyJn(y)^_1(x)-nh2  (x)]   (18) 
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where  y  is  the  absolute  value  of  the  complex  index  of  refraction 

outside  the  sphere,  e  is  the  absolute  value  of  the  complex 

index  of  refraction  inside  the  sphere,  x  =  2frr/A  and  y  =  mx. 

The  functions  i     ,  y  ,  h  2  =  j  -iy   are  the  spherical  Bessel  and 
Jn'-'n    n    JnJn  r 

Neuman  functions.   A  procedure  for  calculating  their  values  is 
outlined  in  the  following  sections. 
3.   Computation  of  j  ,  y  ,  h2 

The  argument  of  these  Bessel  related  functions  is 
complex  in  some  cases,  so  much  care  is  required  in  their  calcu- 
lation.  Each  function  is  calculated  using  recursion  relations 

and  the  lowest  order  functions  as  follows.   For  j  , 

Jn 

(2n+l)j  (z) 

j  .,(z)  =  =-£ -  j   ,(z)  (19) 

Jn+1  Z  Jn-1 

.  ,  v    sinz 

j  (2)  =  sinz  _  cosz  (20) 


j0(z)  =  (  —  -  — )sinz  -  —  cosz 
2         z3   z 


likewise  for  y  , 


n 


(z)  =  ^2n+1>7n^)  _  y   (z)  (21) 

n+1  z         ■'n-l 


,     s    -cosz 

y  Cz)  =  — — - 

J  o  z 


-cosz    sinz 


Yl(z)  =   — j —  C22) 


-3     1  3 

y„(z)  =  (  —  +  — )  cosz  -  —   sinz 

2         z3    z  z2 
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and  these  two  yeild  h  z, 

n 


h  2(z)  =  j  (z)-iy  (z)  (23) 

n       Jn     J  n 


Because  the  argument  of  the  sine  and  cosine  in  evaluating  the 
spherical  Bessel  function  is  complex,  the  following  relations 
are  needed, 

z  =  a  +  ib 

cosz  =  cos (a)cosh(b)  -  isin(a)sinh(b)     (2H) 

sinz  =  sin(a)cosh(b)  +  icos (a) sinh(b) 

where  cosh  and  sinh  are  the  usual  hyperbolic  trigonometric 
functions . 

D.   DESCRIPTION  AND  DOCUMENTATION  OF  PROGRAM 

The  program  adapting  Mie  theory  to  machine  computation  is 
composed  of  the  MAIN  routine  and  numerous  subroutines  described 
in  the  following  sub-sections. 

1.   MAIN  Routine 

The  MAIN  routine  handles  the  input  and  output  functions 
necessary  in  program  operation.   The  input  parameters  include 
indices  of  refraction  inside  and  outside  of  the  spheres,  con- 
stants of  the  particle  distribution  and  the  most  frequently 
occurring  value  of  Mie  size  parameter  in  the  distribution,  x  . 
Through  x  ,  the  wavelength  of  the  incident  light  is  entered 
because  r   is  known  for  any  desired  particle  distribution. 
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Other  necessary  inputs  are  the  number  of  scatter  angles  desired 
in  the  output  scattering  matrix  elements  P.(0),  the  smallest 
value  of  x,  the  increment  in  x  to  be  used  for  numerical  inte- 
gration and  the  odd  number  of  x  values  at  which  the  integrand 
is  to  be  evaluated.   The  MAIN  program  accepts  the  input  values 
and  uses  the  distribution  parameters  to  assign  each  designated 
x  value  a  weight.   Because  the  particle  distribution  is  normali- 
zed, the  sum  of  all  weights  assigned  is  unity. 

Simpson's  1/3  rule  is  used  to  evaluate  integrals  5,  6 
and  7.   At  each  value  of  x  the  scattering  and  extinction 
efficiency  factors  of  equations  11  and  12  are  calculated  through 
use  of  a  subroutine  MIEM.   Similarly,  at  each  desired  scatter- 
ing angle  the  dimensionless  intensity  parameters  of  13  are 
calculated  through  MIEM.   Printout  can  be  called  at  each  value 
of  x  if  so  desired.   The  MAIN  program  cumulatively  sums  the  areas 
of  an  even  number  of  panels  to  get  the  desired  results  which  are 
then  printed. 

2.   MIEM  Subroutine 

MIEM  is  a  subroutine  composed  to  calculate  the  scatter- 
ing and  extinction  efficiency  factors  in  the  integrand  of 
equations  5  and  6  and  the  dimensionless  intensity  parameters  of 
equation  13.   Required  inputs  of  MIEM,  transferred  to  it  by 
MAIN,  are  the  index  of  refraction  of  the  particle  (henceforth 
the  index  of  refraction  of  the  outside  medium  will  be  unity) , 
the  size  parameter  x,  and  the  number  of  scattering  angles  at 
which  calculation  of  the  intensity  parameters  is  desired.   MIEM 
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uses  equations  11  and  12  to  calculate  the  efficiency  factors 

and  equation  13  to  calculate  the  intensity  parameters. 

Because  equations  9,  10,  11  and  12  require  summation  of  series, 

each  element  of  the  series  must  be  evaluated  in  turn.   The 

terms  involved  in  finding  the  efficiency  factors  are  a   and 

&  3  n 

b   which  are  calculated  for  each  n  =  1,  2,  3,  4, until 

n  '   '   ' 

the  next  term  of  the  series  is  adequately  small  or  the  total 

number  of  terms  exceeds  120.   After  each  calculation  of  a   and 

b   its  contribution  to  Q    (x)  and  Q   , (x)  is  added  to  the  pre- 
n  xsca         ext 

vious  total.   MIEM  uses  function  subroutines  JN ,  FN  and  HN  to 

calculate  each  of  the  a  's  and  b  's.   Arrays  JX,  FX  and  HX  are 

n        n         j  •> 

used  to  store  values  of  the  corresponding  spherical  Bessel 

functions  during  each  recursion  step. 

After  all  the  parameters  a   and  b   are  calculated, 

n      n 

MIEM  turns  them  over  to  ANGLE  to  complete  the  remaining  dimension- 
less  intensity  parameter  calculations. 
3.   ANGLE  Subroutine 

ANGLE  is  a  subroutine  called  by  MIEM  to  calculate  the 
dimensionless  intensity  parameters  at  each  desired  angle  of 
scatter.   ANGLE  requires  as  inputs  the  Mie  scattering  coeffici- 
ents generated  in  MIEM  and  the  number  of  scatter  angles  desired. 
When  called,  ANGLE  uses  equations  9  and  10  to  find  S,  and  S_  at 

each  scatter  angle.   To  do  so  ANGLE  evaluates  t   and  tt   using 

°  n      n      ° 

recursion  relations  14  and  15  with  initial  order  functions,  16. 

At  each  n,  the  corresponding  Mie  coefficients  a   and  b   are  used 
'  x-      o  n      n 

with  t   and  tt   and  their  contributions  to  Sn  and  S0  are  added, 
n      n  12 
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The  process  is  repeated  for  each  scatter  angle  and  the  result- 
ing arrays  are  used  to  evaluate  equation  13  at  each  scatter 
angle.   The  resulting  dimensionless  intensity  parameters  are 
then  returned  to  MIEM  and  in  turn  to  MAIN  for  integration  using 
equation  7 . 

4.   JN,  FN,  HN,  CCOS,  CSIN  Complex  Functions 

These  functions  are  called  by  MIEM  in  evaluation  of 

the  Mie  coefficients  a   and  b  .   JN ,  FN  and  HN  contain  logic 

n      n      '  & 

to  perform  the  recursion  operation  of  equations  19  through  23. 
CCOS  and  CSIN  are  complex  trigonometric  functions  drawn  upon 
as  needed  by  JN,  FN  and  HN. 

E.   CPU  TIME  CONSIDERATIONS 

CPU  time  depends  largely  on  the  wavelength  to  particle  size 
ratio.   This,  of  course,  varies  with  the  particle  size  distribu- 
tion used.   Thus,  if  the  distribution  includes  many  particles 
of  large  size  compared  to  the  wavelength  the  CPU  time  is  great 
and  vice  versa.   Typical  CPU  time  requirements  were  on  the  order 
of  15  to  30  minutes  for  wavelengths  of  .53  to  .28  microns,  with 
about  300  x  values  in  a  distribution  of  Water  Cloud  C.2  of 
Deirmendjian  [21].   The  dependence  on  particle  size  to  wave- 
length ratio  is  due  to  the  fact  that  many  terms  are  required  for 
convergence  of  series  for  large  values  of  x.   As  expected,  CPU 
time  goes  up  quite  linearly  with  the  number  of  x  values  used  as 
integration  points.   There  is  little  dependence  on  the  number  of 
scatter  angles  required. 
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APPENDIX  C 

I,   SIMULATION  OF  A  CLOUD  IN 
MONTE  CARLO  ROUTINE  LITE 

A.   INTRODUCTION 

Reference  6  describes  Monte  Carlo  simulation  of  light  propa- 
gation through  an  infinite,  homogeneous  atmosphere.   Many 
problems  can  be  sufficently  investigated  using  this  model,  but 
one  of  the  advantages  of  a  Monte  Carlo  model  is  its  relative 
ease  in  adaptation  to  irregular  geometries  and  inhomogeneous 
atmospheres.   This  Appendix  gives  one  method  for  simulation  of 
a  cloud  with  plane  parallel  homogeneous  medium.   All  macro- 
scopic properties  are  the  same  everywhere  inside  of  the  cloud 
and  another  set  of  macroscopic  parameters  are  the  same  every- 
where outside  of  the  cloud.   This  model  circumvents  the  entire 
problem  of  describing  and  locating  boundaries  and  inhomogeneities 
in  real  clouds.   Figure  5   depicts  the  general  structure  of  the 
cloud  model  giving  the  names  of  various  parameters  of  the  com- 
puter simulation.   A  point  source  of  light  is  incident  normal 
to  the  left  edge  of  the  cloud  and  its  path  is  randomly  generated 
until  it  exits  the  outermost  sphere  of  the  reference  volume . 
Reference  6  describes  the  random  path  generation  and  accounta- 
bility also  used  in  this  model  so  the  terminology  of  that 
reference  is  used  here  for  continuity  whenever  possible. 
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B.  MODELING  CLOUD  BOUNDARIES  AND  PROPERTIES 

The  sets  Qf  parameters  needed  in  defining  two  different 
scattering  media  are  implemented  by  use  of  storage  arrays .   A 
binary  logic  code  is  used  that  switches  whenever  a  photon 
crosses  a  boundary.   Each  array  of  characteristic  parameters 
has  two  columns,  one  for  inside  the  cloud  and  one  for  outside 
the  cloud.   The  logic  switch  determines  at  each  collision 
from  which  column  the  parameter  is  to  be  drawn. 

Upon  crossing  a  boundary,  the  photon  is  stopped  and  a  new 
distance  is  randomly  generated  in  accordance  with  the  para- 
meters of  the  newly  entered  medium.   The  method  used  for 
determination  of  whether  or  not  a  boundary  was  crossed  is 
described  in  the  following  section. 

C.  DETERMINATION  OF  BOUNDARY  CROSSING 

The  following  equation  is  used  to  determine  at  each  colli- 
sion the  angle  between  the  original  direction  of  incidence  and 
the  present  position  vector,  R, 

-i  *-(*-v 

0  =  COS      


where  r  is  the  distance  between  the  collision  point  and  the 
point  of  incidence  on  the  cloud.   R,  is  an  orientation  vector 
from  the  point  of  collision  to  the  CO,  0,  1)  point  of  the  fixed 
coordinate  system.   The  r  and  0  values  are  computed  at  each 
collision  which  allows,  due  to  cylindrical  symmetry,  calculation 
of  the  projected  distance  along  the  direction  of  incidence. 
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This  distance  is 

Projected  distance  =  r  cos  0  =  DIST 

new 

which  is  compared  to  the  thickness  of  the  cloud.   The  medium 
in  which  the  photon  exists  is  determined  by, 

DIST     <  0  Behind 

new  I 

0  <      DIST     <  THICKNESS) Inside 
new  J 

DIST     >  THICKNESSl Beyond 
new  \  J 

The  location  of  the  collision  relative  to  the  cloud  is  compared 
to  the  location  of  the  previous  collision  and  a  boundary 
crossing  is  found  if  it  has  occurred  between  collisions.   Logic 
was  created  to  then  stop  the  photon  at  the  boundary  and  project 
it  along  the  same  path  using  new  scattering  parameters.   There 
are  nine  different  combinations  of  old  and  new  positions ,  three 
possibilities  for  the  old  position  and  three  possibilities  for 
the  new  position.   Each  specific  situation  is  investigated  at 
each  collision  and  upon  boundary  crossing  the  correct  stopping 
formula  is  applied. 
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APPENDIX  D 

I.   AUTOMATION  OF  RELATIVE  FLUX  CONTOUR 
PLOTTING  IN  DRLITE  ROUTINE 

A.  INTRODUCTION 

One  very  important  output  of  the  Monte  Carlo  routine  is  the 
relative  flux  per  unit  area  at  numerous  grid  points  within  the 
reference  spheres  with  respect  to  the  original  energy  projected 
by  the  point  source.   Knowing  the  relative  flux  at  numerous 
grid  points  allows  generation  of  equal  photon  flux  contours. 
Calculation  of  the  contours  was  a  tedious  job  requiring  many 
hours  on  a  programmable  pocket  calculator  and  at  the  drawing 
board.   Accuracy  of  the  contours  suffered  (not  to  mention  the 
one  doing  the  work)  therefore  a  method  was  created  to  produce 
many  of  the  contour  plots  found  within  this  report. 

B.  AUTOMATION 

A  subroutine,  CONISD,  was  found  in  the  computer  program 
library  that  was  designed  to  produce  a  contour  map  of  irregu- 
larly spread  data  points.   Each  data  point  is  a  triad  of  x,  y 
and  z  values  where 

z  =  f(x,y) 

This  procedure  is  nicely  suited  for  plotting  the  flux  contours 
required  since  the  flux  per  unit  area  has  been  calculated  at 
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each  increment  in  range  and  scattering  angle.   In  transforming 
evenly  spaced  range  and  angle  intervals  from  polar  to  cartesian 
coordinates,  irregularly  spaced  intervals  are  formed.   By 
inspection  the  flux  varies  quite  linearly  in  the  range  direction, 

The  region  to  be  contoured  is  subdivided  into  triangular 
areas,  using  acute  triangles  as  much  as  possible.   The  method 
of  triangularization  is  based  on  that  of  Ref.  30.   The  list  of 
triangles  is  then  analyzed  for  adjacencies.   For  a  given  con- 
tour level,  each  outer  boundary  of  the  area  to  be  contoured 
is  checked  for  a  possible  entry  value.   Upon  finding  a  value, 
the  contour  line  is  traced  from  triangle  to  triangle  until  it 
exits  the  area.   As  the  line  is  traced  from  segment  to  segment, 
the  four  nearest  values  are  linearly  interpolated  and  the 
resulting  value  stored.   After  all  lines  have  been  found  and 
intersections  stored,  another  subroutine  is  called  to  fit 
smooth  curves  to  the  stored  values.   The  resulting  curves  are 
the  desired  contours  of  equal  photon  flux. 

C.   ADDITIONS  TO  DRLITE 

A  section  was  created  within  DRLITE  to  drive  the  subroutine 
CONISD  after  computation  of  relative  flux  at  desired  grid 
points.   Inputs  required  by  CONISD  are  the  number  of  grid  points, 
each  corresponding  triad  of  values ,  the  value  of  each  contour 
desired,  scaling  factors  for  plotting  the  output  and  designa- 
tion of  which  data  points  are  to  be  boundary  points . 

Because  data  used  to  generate  the  contours  was  generated  by 
a  Monte  Carlo  routine,  a  measure  of  statistical  significance  had 
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to  be  included  before  actual  plotting.   The  measure  used  was 
based  purely  on  the  number  of  photon  crossings  at  each  parti- 
cular angular  bin  at  every  range.   If  the  number  of  crossings 
was  less  than  three,  the  corresponding  data  triad  was  removed 
from  the  list  to  be  used  in  generation  of  contours. 

D.   DISCUSSION 

The  automation  method  explained  above,  although  useful  in 
reducing  manual  labor,  is  quite  consumptive  of  CPU  time.   Com- 
puter time  goes  up  almost  expotentially  with  the  number  of 
grip  points  used  in  the  routine.   This  can  be  attributed  to 
the  number  of  searches  necessary  in  passing  from  segment  to 
segment.   Although  the  general  idea  behind  automation  of  the 
'contour  plotting  is  good,  the  actual  implementation  of  such 
methods  must  be  done  with  caution. 

As  an  aside,  the  subroutine  CONISD  also  has  the  capability 
to  have  internal  cut  out  boundaries  placed  into  its  logic. 
Future  application  of  this  feature  may  be  to  draw  equal  photon 
flux  contours  for  light  propagation  around  defined  objects. 
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APPENDIX  E 

I.   DERIVATION,  DOCUMENTATION  AND  VERIFICATION 
OF  EFFECTIVE  ATTENUATION  COEFFICIENT  METHOD 


A.  INTRODUCTION 

Gordon  [9]  presents  the  concept  of  an  effective  attenuation 
coefficient  (EAC)  which  considerably  simplifies  multiple  scat- 
tering calculations.   Concise  expressions  for  the  total  flux 
through  an  aperture  and  the  beam-spread  function  are  derived 
in  terms  of  the  EAC.   A  program  was  written  to  evaluate  these 
expressions  so  that  comparison  with  other  light  scattering 
models  could  be  accomplished.   This  Appendix  derives  the 
expressions  given  by  Gordon,  documents  the  routine  written  to 
evaluate  the  expressions  and  shows  verification  of  program 
accuracy. 

B.  DERIVATION 

Gordon  first  suspected  that  simple  formulas  might  describe 
multiple  scattering  after  examining  Duntley's  [29]  measurements 
of  the  fraction  of  power  emitted  from  a  collimated  source  which 
reaches  a  circular  aperture  subtending  a  cone  of  half-angle  9 
when  viewed  from  the  source  and  located  a  distance  r  from  it. 
Gordon  noticed  that  on  semi-log  paper  the  data  formed  straight 
lines  to  about  15  extinction  lengths  when  plotted  as  a  function 
of  range.   The  expression  for  flux  reaching  the  aperture  where 
F   is  the  source  strength  can  be  expressed  as  [9]: 
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F(0,R)  a  F  e 


-ap(9)R 


CI) 


where  a e(8)  is  the  EAC  as  a  function  of  0. 

Consider  the  geometry  of  Figure  38.   A  unidirectional  point 
source  is  centrally  located  on  the  plane  x  =  0.   On  the  plane 
x  =  R  is  located  a  circular  aperture  which  subtends  a  cone  of 
half-angle  9  with  the  point  source.   The  total  flux  through 
the  aperture  is  given  by  [9]: 


where , 


F  =  F  exp 
o  * 


J  ,     2ti  9'  v 

-11--"-/  I  I  I    P(9)dJ  dtJoR 


•V    *  0  0 


(2) 


doj  =  sin  9  d9  d0  =  differential  solid  angle 
s  =  scattering  coefficient  in  km 


«  =  extinction  coefficient  in  kn 
t  =  x/R  =  normalized  range 


-1 


0  *  =  tan 


•l/tan  Q\ 
ll-x/RJ 


and  P(9)  is  the  normalized  phase  function,  thus  it  satisfies  the 
relation, 


ff 

oo 


P(9)  dco  =  1 


(3) 
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Figure  3  8 
Geometry  of  EAC  Method 
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A  simple  physical  description  of  equation  2  is  that  it 
accounts  for  flux  absorbed  and  scattered  from  the  beam  in  the 
first  term  of  the  exponent  and  also  accounts  for  flux  that  is 
scattered  back  through  the  aperture  in  the  integral  term  of 
the  exponent . 

Comparing  equation  1  with  equation  2  the  effective  attenu- 
ation coefficient  is  found  to  be 

1   2tt  9" 
ae(0)  =  a(l-M  J  J  J   PCe)do)|dt]  (4) 

In  adapting  this  equation  to  machine  computation,  a  normali- 
zed phase  function  is  necessary.   A  reasonable  choice  is  the 
Henyey-Greenstein  phase  function  [73, 

P(8)  =  (1~G2) -1<G<1        (5) 

t+Tr(l+G2-2Gcos0)1  '5 

Because  the  Monte  Carlo  model  makes  use  of  this  phase  function 
also,  comparison  of  the  two  models  is  possible. 

Now  that  the  flux  through  an  aperture  is  known,  it  can  be 
used  to  derive  the  irradiance  caused  by  a  unit  strength  uni- 
directional point  source  at  each  point  of  a  selected  target 
plane.   This  irradiance  distribution  is  called  the  beam  spread 
function,  HR(0,R).   Let  the  initial  point  source  have  unit 
strength.   By  considering  the  flux  through  a  small  annulus  of 
width  R  d0  located  at  (R,0),  HD(0,R)  can  be  written  as  [9], 
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HB(6,R) 


1  dF 

2TrR^sine         d6 


(6) 


where  dF/d6  is  the  derivative  of  equation  2  with  F   =1.   To 
evaluate  this  derivative,  Leibnitz  rule  is  used  [33].   It  states 
that  if 


=  /  0(x, 
J(+.) 


F(t)  =|   0(x,t)ctt 
(t) 


(7) 


where  a(t)  and  b(t)  are  dif ferentiable  functions  of  t  and  where 

3  0 
0(x,t)  and  ~ —  are  continuous  in  x  and  t,  then 


3t 


(t) 


dF 
dt 


/80(x,t 
(t)  3t 


I   dx+0[b(t),t]d^t:>  -  0[a(t),t3d^lt)  (8) 


dt 


dt 


Take  the  derivative  of  equation  2 


dF    d_ 
d9  "   d6 


1   2tt 


-i-l 


(tffr^h) 


2"FsR  3F1 


/' 


P(6)sin  9  d0 


dt 


(9) 


(10) 


and  apply  Leibnitz  rule  where 


,e)  =J 


0(t,6)  =  /   P(9)sin  0  d8 

-'o 


a(t) 
b(t) 


0 
1 


(11) 


it  is  now  apparent  that  the  derivative  becomes, 
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§-  2^fsr/(tI/p(9)  sin 


0  d0  dt 


(12) 


The  derivative  that  remains  is  of  an  integral  whose  range  is 
over  the  same  independent  variable,  6,  therefore  the  result  is 
the  product  of  the  integrand  evaluated  at  the  upper  limit  and 
the  derivative  of  the  upper  limit,  i.e., 


r/« 


!^/p(9>  sin6d9  =  P(0')  sin(6')  |^-(ef)    (13) 


but 


3_ 
89 


tan 


-1  /tan( 


1-t 


i+pey 


(14) 


(l-t)cos20 


so  it  is  clear  that 


f  =  2*FsR/ 
•'o 


P(9'  )sin(0'  )  -= 


(l-t)cos2! 
Using  the  small  angle  approximation,  tan0  =  0  =  sin0  and  letting 


^tsney 


dt      (15) 


u    =    tan 


■lW 


(16) 


du    =  -= 


dt 


1  + 


1-t 


(17) 


(1-t) 


equation  15  takes  the  form 
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#fir/2 
II     P( 


||  =    2irFrjH  I      P(u)    cos(u)    du  (  I  ••■•  * 


Using  equation  18  in  equation  6  the  beam  spread  function  is 

„ir/2 


hi     ^ 


HbO,R)  =  —-   I       P(u)  cosu  du  (  !  :i  ) 


A  useful  calculation  in  comparing  this  beam  spread  function 
with  that  found  using  the  Monte  Carlo  method  is  taking  the 
relative  logarithm  of   the  beam  spread  function. 

C.  ADAPTATION  OF  EFFECTIVE  ATTENTUATION  COEFFICIENT  METHOD  TO 

MACHINE  COMPUTATION 

The  Effective  Attentuation  Coefficient  method  lends  itself 
nicely  to  machine  computation.   The  inputs  required  by  the 
program  are  the  scattering  coefficient,  the  extinction  coef- 
ficient and  the  Henyey-Greenstein  G  parameter.   Also  entered 
are  the  initial  range  and  theta,  the  increment  in  range  and 
theta  and  the  number  of  grid  points  desired  in  each  dimension. 

The  routine  written  is  called  EAC .   EAC  calls  on  a  library 
subroutine  DQG32  which  integrates  functions  FCC  and  FCT  inside 
of  which  the  integrands  of  equation  4-  and  19  are  defined, 
respectively.   DQG32  uses  32  point  Gaussian  quadrature  which 
integrates  polynomials  up  to  degree  6  3  exactly. 

The  output  of  EAC  can  be  received  by  numerical  matrix  tabu- 
lation or  by  graphical  contour  plotting  as  described  in  Appendix 

D.  Examples  of  these  plots  are  found  in  Figures  6  through  9. 
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D.   VERIFICATION  OF  EAC 

The  results  from  EAC,  the  NPS  routine,  were  compared  to 
figures  obtained  by  Gordon  [9]  using  the  same  method.   There  is 
a  small  dependence  on  phase  function  for  small  optical  thick- 
nesses for  which  this  method  is  applicable.   Thus  in  comparing 
results  with  Gordon,  a  slight  deviation  is  noticed  due  to  the 
fact  that  the  phase  function  used  may  not  have  been  exactly  the 
same.   A  Henyey-Greenstein  G  parameter  of  .94  was  used  to  gen- 
erate the  figures  presented  here.   Figure  39  compares  the 
results  of  EAC  with  Gordon  for  irradiance  as  a  function  of  6 
at  a  range  of  2.12  extinction  lengths.   Agreement  is  seen  to  be 
quite  good.   Figure  40  compares  the  results  of  EAC  with  Gordon 
for  irradiance  as  a  function  of  range  with  a  scatter  angle  of 
10  degrees.   Agreement  is  good  here  also.   These  are  a  few  of 
many  checks  on  routine  accuracy  which  lead  to  much  confidence 
in  the  routine. 
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Verification  of  EAC  Method 
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APPENDIX  F 
I.   CHECKS  ON  POSSIBLE  ERRORS 

A.  GENERAL 
Each  computational  procedure  employed  in  this  simulation 

was  verified  when  possible  and  every  effort  was  made  to  ensure 
the  entire  simulation  functioned  properly.   Results  of  various 
theories  were  compared  using  identical  parameters  to  build 
confidence  at  intermediate  steps  along  the  path  toward  final 
results.   Procedures  used  to  document  the  correct  behavior 
of  the  computer  routines  used  in  this  simulation  are  described 
in  this  Appendix. 

B.  SPECIFIC  CHECKS 

One  of  the  secondary  objectives  of  this  report  was  to  com- 
pare results  of  different  theories  in  predicting  light  transfer 
through  a  scattering  absorbing  medium.   The  ability  to  compare 
results  of  separate  theories  stands  by  itself  as  a  verification 
of  the  accuracy  of  each  theory. 

Reference  6  states  numerous  initial  methods  used  to  verify 
the  mechanics  of  the  computer  simulation  in  its  first  stages. 
In  converting  the  routine  from  one  that  simulates  an  infinite 
homogeneous  medium  to  one  that  contains  a  cloud,  additional 
checks  were  required.   Many  photons  were  traced  manually  by 
hand  calculations  to  ensure  that  all  nine  possible  combinations 
of  transfer  across  boundaries  were  correctly  handled  by  the 
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routine.   Errors  were  found  and  corrected  on  many  an  occasion. 
In  the  situation  where  the  transmission  parameters  and  the 
phase  function  of  the  cloud  were  the  same  as  those  of  the 
surrounding  medium,  one  would  expect  the  results  to  be  identi- 
cal to  those  found  when  simulating  a  homogeneous  medium.   This 
in  fact  was  the  case.   Results  were  identical  to  well  within 
statistical  scatter. 

Changeover  from  measurements  in  terms  of  extinction  lengths 
to  measurements  in  terms  of  real  dimensions  also  required  brief 
checking.   One  would  expect  the  results  to  be  identical  when 
comparing  a  medium  with  an  extinction  coefficient  of  one  in- 
verse kilometer  to  a  medium  described  in  terms  of  optical 
thickness.   This  is  so  because  one  unit  of  distance  is  the 
same  in  both  cases,  one  kilometer.   This  is  the  case.   In  com- 
paring trials  where  the  first  has  a  distance  between  accounta- 
bility shells  of  twice  the  second  trial  but  an  extinction 
coefficient  of  half  the  second  trial  with  the  albedo  the  same, 
the  number  of  photons  crossing  each  bin  should  be  the  same. 
Once  again  this  is  the  case. 

In  some  cases  the  results  were  checked  by  intuition  only 
as  no  previous  work  with  similar  presentation  was  found. 
Specifically,  the  equal  photon  flux  contours  when  passing 
through  clouds  were  checked  only  by  what  one  would  expect  them 
to  look  like. 

Many  test  runs  were  made  using  the  model  generating  a 
scatter  angle  weighted  by  an  arbitrary  phase  function.   The 
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angles  generated  are  in  fact  weighted  by  areas  of  panels  and 
linearly  within,  panels. 

Confidence  in  the  statistical  nature  of  a  Monte  Carlo 
routine  is  related  to  the  number  of  photon  histories  used  in 
each  trial.  The  smallest  number  of  photons  possible  were  used 
to  create  outputs  of  the  desired  accuracy  in  order  to  minimize 
computer  time  consumption.   In  most  cases  of  spatial  spread 
thousands  of  photon  histories  were  generated.   In  time  spread 
cases  as  few  as  100  were  required  at  times. 

Other  checks  made  are  found  in  Ref.  6  as  well  as  within 
the  body  of  this  work. 
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APPENDIX  G 

I.   DERIVATION  OF  CLOSED  FORM  EXPRESSION 
FOR  TIME  SPREAD 

Stotts  [26]  derives  a  closed  form  expression  for  time  spread 
as  follows .   From  vector  analysis ,  the  average  pathlength  of  a 
photon  per  unit  time  is 


dR 
dt 


(^)2  +  (dj)2  +  (dZ)2 
dt     dt     dt 


which  reduces  easily,  where  r  =  xl  =  yj ,  to 


dR 
dt 


1  + 


dr  I  2 
dz1 


dz 


From  basic  scattering  theory  [31]  the  projected  angle  of  scatter- 
ing is  approximated  by 


[*«*  Y^ 


where  wQ is  the  albedo  of  single  scattering,     t   is  the  optical 
thickness  of  the  scattering  region  and  yQ   is  the  RMS  scatter 
angle.   Thus  the  projected  mean-squared  transverse  displacement 
is 


dr  |  =  z(0)ox  Yq)  2 
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thus 


and 


1 dz1    4   o   'o 


dR  =  [1+i  (oqt  Yq]%  dz 


Integrating  results  in 


R  = 


8z 


2  7cj  T  y^ 
o   '  o 


fl+|  0)  x  Y2V5-1 

V  4  o   oy 


The  average  multipath  time  spread  is  defined  as  the  difference 
between  the  average  transient  time  incurred  from  multiple 
scattering  and  the  normal  transient  time  in  the  absence  of 
scattering.   Thus,  Stotts  has  for  time  spread, 


L  =  Z 
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